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INTRODUCTION 

This  report  summarizes  the  research  progress  made  under  contract 
number  F49620-80-C-0032  entitled  "Research  and  Development  in  Speckle 
Imaging".  Analytic  Information  Processing,  Inc.  pursued  this  effort 
for  the  Air  Force  Office  of  Scientific  Research  (AFSC),  United  States 
*  Air  Force.  The  AFOSR  programs  manager  is  Dr.  H.  Radoski  and  the 

principal  investigator  is  Dr.  J.  W.  Sherman. 
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RESEARCH  OBJECTIVES 


The  objectives  of  the  effort  were  to: 

1.  Develop  a  better  understanding  of  the  random  processes  which 
degrade  the  coherence  of  light  propagating  through  the 
atmosphere. 

2.  Develop  a  first  cut  capability  to  process  atmospherically 
degraded  images  to  obtain  diffraction  limited  resolution. 


To  accomplish  these  research  objects,  the  statment  of  work  for  this 
effort  delineated  six  tasks. 

1.  Develop  mathematical  models  for  the  random  speckle  image 
produced  by  atmospheric  turbulence  and  by  recording  effects. 

2.  Develop  equations  for  the  second  order  statistics  of  a  set  of 
speckle  images. 

3.  Evaluate  reconstruction  performance  in  terms  of  algorithm 
complexity  and  reconstruction  accuracy. 

4.  Develop  computational  techniques  from  the  reconstruction 
equations  with  constraints  describing  the  object's 
characteristics. 

5.  Publish  the  results  of  analyses  in  relevant  scientific 
journals  and  at  technical  conferences. 


6.  Coordinate  the  speckle  imaging  techniques  developed  with 
those  of  other  researchers. 


r 


Although  the  feasibility  of  Speckle  Imaging  had  been  shown  previously, 
P  the  models  and  algorithms  used  were  the  simplest  possible.  The 

development  of  improved  models  and  associated  algorithms  was  necessary 
as  a  basis  for  improved  Speckle  Imaging  and  for  analysis  of  the 
accuracy  of  the  reconstructed  images.  These  improved  models  were 
l  developed  by  analytic  and  experimental  means.  These  tasks  emphasized 

the  development  of  models  which  allow  for  the  natural  variation  of  the 
degradation  and  the  development  of  speckle  imaging  algorithms  which 
estimate  not  only  the  underlying  image  but  also  the  model  parameters 
r  or  degradation  conditions. 


» 


> 


» 
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STATUS  OF  RESEARCH  TASKS 


The  status  of  the  research  tasks  will  be  discussed  in  terms  of  two 
areas: 


1.  Non-isoplanatic  Modeling 

2.  Recursive  Algorithms 


We  have  chosen  to  concentrate  our  research  in  these  two  areas  because 
they  appear  to  have  the  greatest  potential  for  improving  the  range  of 
applications  for  Speckle  Imaging.  The  treatment  of  the 
non-isoplanatic  case  will  allow  Speckle  Imaging  to  be  used  over  wider 
fields-of-view  in  order  to  image  extended  objects.  The  development  of 
recursive  algorithms  will  allow  the  use  of  generally  available 
mini -computers  to  process  Speckle  Images  because  they  reduce  the  large 
I/O  requirements  of  non-recursive  algorithms. 

Non-Isoplanatic  Modeling 

Our  modeling  efforts  for  Speckle  Imaging  under  non-isoplanatic 
conditions  were  designed  to  determine  the  nature  of  the  degradation  of 
Speckle  Images  and  the  extent  to  which  the  underlying  image  can  be 
recovered.  The  two-point  source  wave-structure  function  and  the 
altitude  dependence  of  the  index-of-refraction  structure  constant  were 
used  as  the  basis  of  the  modeling.  The  first  and  second  order 
statistics  of  speckle  images  were  derived.  An  analysis  of  the 
information  recoverable  from  these  statistics  has  shown  that 
diffraction  limited  resolution  can  be  obtained.  Two  important 
assumptions,  which  are  usually  satisfied,  were  made: 

1.  The  Isoplanatic  patch  encompasses  the  average  spread 
function. 
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2.  The  greatest  contribution  to  the  atmospheric  distortion  comes 
from  turbulence  at  low  altitudes. 

Under  these  conditions,  the  aberrations  of  the  telescope  do  not  affect 
the  second-order  statistics.  Also,  the  integral  equation  relating  the 
object  and  measured  second-order  statistics  is  well  conditioned. 

These  results  were  presented  at  the  SPIE  annual  symposium  and 
published  in  the  SPIE  Proceedings,  Volume  243. 

This  modeling  effort  was  continued.  The  parameters  of  the 
mathematical  model  for  the  non-isoplanatic  degradation  are  determined 
from  the  i ndex-of-refraction  structure  constant  profile.  A  random 
process  model  of  the  index-of-refraction  structure  constant  profile 
was  used  to  derive  the  dependence  on  this  profile  of  the  function 
which  smooths  the  second-order  statistics  under  non-isoplanatic 
conditions.  We  found  that  for  the  purposes  of  this  model  the 
atmospheric  structure  function  was  reasonably  approximated  by  a 
5/3-law  model  for  all  values  of  its  argument.  Simple  useful  models 
were  obtained  for  the  smoothing  function  which  describes  the 
non-isoplanatic  conditions.  These  models  were  verified  numerically. 

We  plan  to  submit  these  results  to  the  Journal  of  the  Optical  Society 
of  America  for  publication. 

The  results  of  this  modelling  effort  were  used  to  develop  a 
restoration  algorithm  for  non-isoplanatic  Speckle  Imaging.  The 
tradeoffs  in  computational  complexity  were  examined.  The  Speckle 
Imaging  equations  were  formulated  in  a  general  framework  based  on 
linear  algebra  to  allow  an  examination  of  the  many  alternatives.  We 
plan  to  submit  these  results  to  the  Journal  of  the  Optical  Society  of 
America  and/or  the  IEEE  Transaction  on  ASSP  for  publication. 

Recursive  Algorithms 

The  recursive  algorithm  development  effort  was  pursued  because  it 
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allows  Speckle  Imaging  to  be  used  for  a  wider  range  of  applications. 
Extended  Kalman  filter  formulations  were  developed  so  that  subject 
image  models  can  be  introduced  easily.  For  example,  time  evolutionary 
models  for  solar  features  or  restrictive  models  like  limited  spatial 
extent  can  be  used  easily.  The  recursive  nature  of  the  algorithm  will 
allow  the  use  of  mini -computers  for  data  reduction,  because  of  the 
removal  of  the  requirement  for  storing  the  second-order  statistics. 
These  second-order  statistics  require  millions  of  words  of 
intermediate  storage  which  must  be  accessed  many  times  for 
non-recursive  algorithms,  producing  an  I/O  burden  too  great  for  most 
mini -computer  systems. 

Recursive  algorithms  based  on  an  extended  Kalman  filter  approach  have 
been  derived.  A  detailed  analysis  of  the  decomposition  of  the  Speckle 
Imaging  equations  into  manageable  subsets  has  been  performed.  The 
decomposition  has  a  "natural"  format  that  results  in  very  reasonable 
computation  requirements. 

The  evaluation  of  a  one-dimensional  example  was  started.  An 
oversimpl ication  in  the  models  led  to  singularity  problems.  In 
finding  the  problem,  we  also  found  simpler  forms  for  the  Kalman  gain 
equations.  These  simplified  Kalman  equations  were  compared  to  a 
recursive  form  of  the  least  square  algorithm  previously  developed. 

We  plan  to  submit  these  results  to  the  Journal  of  the  Optical  Society 
of  America  and/or  the  IEEE  transaction  on  Acoustics,  Speech,  and 
Signal  Processing. 
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Section  1  Introduction  and  Definition  of  the  Problem 
1-1  Def ini tion  of  the  Problem 

A  solution  to  the  problem  of  reconstructing  telescopic  images  whose 
resolution  has  been  degraded  by  passage  through  the  turbulent 
atmosphere  has  already  been  proposed  with  the  process  called  "Speckle 
Imaging".  Speckle  Imaging  is  the  reconstruction  of  diffraction 
limited  images  from  atmospherically  degraded,  multi-frame  imagery. 

The  process  of  recording  the  images  viewed  through  the  telescope 
introduces  noise  effects  making  the  treatment  of  the  turbulence 
effects  difficult.  A  model  has  been  developed  for  both  the  turbulence 
effects  and  the  noise  effects.  With  these  two  models  and  a  data  set 
composed  of  images  of  the  object,  a  speckle  imaging  algorithm  provides 
an  estimate  of  the  object.  The  algorithm  previously  reported  [1]  is 
based  on  a  least-squares  method. 

The  goal  of  this  effort  was  to  develop  a  new  algorithm  in  speckle 
imaging  in  order  to  increase  the  range  of  applications.  The  new 
algorithm  should  meet  the  three  following  requirements: 

1)  To  be  a  real-time  algorithm,  treating  each  image  when  it  arrives, 

2)  To  allow  the  treatment  of  time  evolutionary  models  for  both  the 
object  and  the  atmospheric  statistical  character! sties, 

3)  To  allow  the  use  of  a  mini-computer. 

One  of  our  main  concerns  for  the  use  of  a  mini -computer  will  be  the 
reduction  of  data  storage.  The  recursive  nature  of  a  real  time 
algorithm  will  aid  in  this  reduction. 

1-2  -  Speckle  Image  Model 

In  the  effort  to  demonstrate  the  feasability  of  the  least-squares 
algorithm,  several  simplifying  assumptions  on  the  speckle  image  model 
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were  used  and  the  new  algorithm  will  be  developed  under  the  same  set 
of  assumptions. 

-  The  imaging  is  assumed  to  be  isoplanatic  which  means  the 
speckle  image  from  a  point  source  is  constant  over  the  field  of 
view.  This  allows  for  the  use  of  Fourier  techniques  to  convert 
convolutions  into  algebraic  equations, 

-  very  simplified  mathematical  models  for  the  statistics  of  the 
point  spread  response  are  assumed, 

-  to  avoid  any  analysis  or  question  of  degrading  effects,  the 
imaging  conditions  are  restricted  to  narrowband  exposures  of  a 
constant,  bright  object. 

Models  for  the  two  principal  sources  of  image  distortion  that  occur 
when  using  a  ground-based  telescope  can  be  developed  under  the 
previous  assumptions.  The  first  source  is  the  turbulence  of  the 
intervening  atmosphere.  The  second  source  is  the  error  or  noise 
associated  with  the  physical  measurement  of  the  image  and  we  shall 
refer  to  it  as  the  noise  of  measurement.  This  latter  affects  the 
quality  of  the  imaging  less  than  the  former,  but  it  is  its  presence 
which  has  previously  made  the  correction  for  the  turbulence  difficult. 

Model  for  the  Turbulence  Effects : 

For  simplicity,  consider  incoherent  imaging  over  an  isoplanatic  patch. 
The  image  can  be  expressed  as  a  convolution  of  the  object  and  a  point 
spread  response  function. 

if  Uyt)  =  fj  0(h)  s(&-& }t)d2  (1) 

_  oO 

where  y  (_x,t)  is  the  formed  image,  o(^)  is  the  object,  and  sU,t)  is 
the  point  spread  response.  The  point  spread  response  depends  on  the 
optical  system  and  the  phase  delay  and  amplitude  disturbance  processes 
caused  by  turbulence.  Let  S(z,t)  be  the  optical  transfer  function 
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corresponding  to  s(x,t),  that  is, 

c*> 

S(ijt)-  Fj  [sfi,  eirfi'ty?'2) i&0>)  (2) 

e.ip[Jlg->l,t)  y  fl/Z  <l,1  Flip,  t)  r)j 

where 

a(  •  )  is  the  aperture  function  of  the  telescope, 

0(  »  )  is  the  phase  aberration  function  of  the  telescope 

Jl'Z  )  is  the  log-amplitude  disturbance  at  the  aperture  caused  by 
turbulence, 

</(*,')  is  the  phase  delay  disturbance  at  the  aperture  caused  by 
turbulence, 

and  «  J  indicates  Fourier  transform  with  respect  to  _x. 

The  random  properties  of  the  function  -SiZ^tJwill  be  characterized  by 
its  mean  S/z)  and  covariance  ^(z.  z2). 

Assuming  the  usual  Gaussian  distribution  for  1  and  / allows  the  mean 
and  covariance  of  S(_z,t)  to  be  written  in  terms  of  the  wave  structure 
function  D(»).  Using  the  results  of  Fante  [2],  the  mean  and 
covariance  of  SU,t)  are  well  approximated  by 

oO 

Sis )  =  evpt  ' Jj<Ll£',Z)A(fi)tTLp[j  (3) 

—  oc> 

and 


o£>  oO 


For  the  mean  S(z)t  the  effects  of  the  atmosphere  and  of  the 
telescope  separate  Into  a  product  form.  SlZ)  is  dominated  by  the 
exponential  factor  and  decays  to  zero  rapidly  as  the  spatial  frequency 
increases  in  magnitude.  For  the  covariance  S5{2tJZ^),  the  expression 
(4)  is  quite  complex. 

However,  the  Fourier  transforms  of  astronomical  data  were  studied  and 
the  analysis  of  their  first  and  second  order  moments  have  shown  that 
the  characteristics  of  the  statistics  of  the  point  spread  response  are 
reasonably  simple  in  spite  of  the  complex  nature  of  (3)  and  (4). 
Therefore,  ad  hoc  approximations  to  (3)  and  (4)  were  developed;  these 
approximations  have  the  same  general  characteristics  as  the  data. 

The  mean  Si? )  was  approximated  by  a  gaussian  function, 


whose  width  was  proportional  to  the  correlation  distance,  rc  ,  of  the 
complex  disturbance  in  the  aperture  plane  resulting  from  the 
turbulence.  rris  small  compared  to  the  diameter  of  the  aperture  D. 

The  approximation  of  the  covariance  was  proportional  to  the  product  of 
two  factors. 

-  one  that  models  the  narrowband  dimensions  of  the  covariance  and 
the  effects  of  decorrelation  across  the  aperture: 


This  term  exponentially  decreases  with  increasing  distance  from 
the  plane 

-  another  one  that  models  the  broadband  dimension  of  the 
covariance: 
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The  -S  (  jL?)  term  models  the  effect  of  removing  the  mean  in 
the  covariance.  The  remainder  of  this  second  factor  approximates 
the  ideal  diffraction  limited  optical  transfer  function. 

Let  the  covariance  of  the  optical  transfer  function  be  denoted  by 

%(2n2>)  =  E  h  S(z_2)fjw 

Because  of  the  symmetry  properties  of  the  Fourier  transforms  of  real 
functions,  it  is  also  true  that: 


Model  for  the  Noise  Effects 


Although  the  formed  images  of  (1)  are  instantaneous 

functions  of  time,  there  is  some  finite  exposure  (i.e.,  averaging) 
time  involved  in  their  measurement.  We  assume  that  these  times  zJTare 
short  enough  to  "freeze"  the  turbulence  effects  and  this  will  allow  us 
to  apply  (1)  for  these  time-average  images: 

-7*47 


and  to  obtain: 


e*= 


r-4i 

;  s 


ystijt 


(8) 


iUTj)--  ff  oil) 


(9) 
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Telescopic  images  are  generally  recorded  using  either  photographic 
film  or  electronic  cameras,  often  at  limited  photon  flux  levels.  For 
both  systems  and  under  various  hypotheses,  the  effect  of  these  noises 
of  measurement  can  be  modeled  using  combinations  of  Poisson,  additive, 
and  multiplicative  noise  terms. 
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Each  type  of  noise  models  physical  random  processes  that  occur  in  the 
measurement  process.  In  the  recording  of  speckle  images  a  particular 
order  is  useful;  additive  noise  modeling  scatted  and  background  light 
occurs  first,  then  Poisson  noise  modeling  photon  reception,  and 
finally  multiplicative  noise  modeling  sensor  non-linearities. 


Then,  if  denotes  the  measured  image,  the  following  diagram  and 

equations  will  represent  the  measurement  noise  process. 


Additive 


Noise 

Poisson 

Process 

Multipl icati  ve 
Noise 


‘faV  =  U0) 


vAV  = 

c 

f 

y  (-7)}  ~  %  7  7  nrr)^' 


The  noise  process  has  a  zero  mean  such  that 


(11) 


Let 


be  an  adcJitive  noise  with  zero  mean  and  variance  of 
!lp(x  t)  be  a  Poisson  process  with  means  np/i) 7  /  and  variance 
0~P  (*)  iJJ,  and 

be  a  multiplicative  noise  with  unity  mean  and  variance  0^  . 

The  mean  and  variance  of  the  Poisson  process  are  given  by 


% 
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where 


r0  is  the  total  photon  flux  detected,  and 

J0  is  the  integral  or  total  expected  image  intensity. 

The  three  measurement  noises  are  all  assumed  to  be  statistically 
independent  among  themselves  and  independent  of  i(l  P~J.  They  are 
assumed  to  be  spatially  and  temporally  uncorrelated. 

1-3  -  The  Restoration  Algorithm 


Sample  Image  Statistics 

For  the  speckle  imaging  algorithm,  the  sample  first  and  second  order 
statistics  of  a  set  of  images  with  measurement  noise  will  be  used  and 
their  relationships  to  the  mean  and  covariance  of  the  point  spread 
response  have  to  be  established. 

First  Order  Sample  Stati sties:  The  set  of  images  is  composed  of  N 
images  taken  at  times,  T  !/■_,*',  then  the  average  of  the  N  images  gives 
the  first  order  sample  statistics: 

—  I  N  \ 

K*)  -  ~7T  ^  (13) 


Because  of  the  zero  mean  of  the  noise  of  measurement,  the  sample  mean 
is  unbiased 


El[C(i  )  J  -  tU) 


and  when  taking  the  expected  value  of  (9),  i(x)  is  given  by: 

oO  _  — 

L(l  )  -  ff  oQ  )  S  J) 


Fourier  transforming  (15)  changes  this  integral  equation  into  an 


algebraic  one 


J (i)~  0(2)  S(2),  (16) 

where  0(1 )  is  the  transform  of  the  object 

xA  j T(z)  is  the  mean  of  the  Fourier  transforms  of  the  images 
If  H(z)  is  the  sample  mean  of  the  Fourier  transform  of  the  N 
images;  then  its  expected  value  is  given  by 

£[!(*)]=  J(i)^  G(z)Scz).  a?) 

Unfortunately,  these  first  order  statistics  contain  very  little  high 
spatial -frequency  information  because  has  a  very  low  pass 

nature.  Because  of  this  low-pass  nature,  the  inversion  of  (17)  is 
ill-conditioned:  small  variations  in  1(2)  would  result  in  large 
variations  in  the  estimation  of  O(l)  and  thus,  the  effect  of  the 
measurement  noise  would  be  increased. 

However,  the  broad-band  nature  of  the  covariance  will  allow 

us  to  use  the  second  order  statistics  to  restore  the  images  in  the 
high  frequency  domain. 

Second  Order  Sample  Stati sties 


The  following  notation  will  be  used  for  the  second  order  statistics: 

-  The  inverse  Fourier  transform  of  the  covariance  of 

the  transfer  function  S{z),  will  give  Kj- ( x./ , x.^) »  the  covariance 
of  the  point  spread  response  s(x.). 

-  The  covariance  of  a  set  of  images  in  the  spatial  domain 
(respectively  in  the  frequency  domain )  is  given  by  R/U, *_x > 
(respectively  S.U  ,z  )  )  and  the  sample  covariance  by  R.(x  ,x) 

/\c  1  1  t  J. 

(respectively  ^  (z^z^)  ). 
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Taking  the  covariance  of  (9)  gives 

Ri(xi}x;L)r  ]LlUl/)-u'xJ)J  J  {18) 

ao 

-  JJ jJoQj  oQj  RsUrhj  *irh)d  Ji  dJi 

-cO 


In  this  integral  equation,  the  object  oQ)  appears  as  the  cartesian 
product  with  itself.  An  unambiguous  solution  of  this  non-linear 
integral  equation  will  require  a  non-linear  solution  method.  The 
method  will  depend  on  the  characteristics  of  the^sample  image 
covariance  used  in  the  inversion  of  (18).  Let  be  the  sample 

covariance  function  of  N  images  at  times  T^.j  =  1,  ...  ,N,  then 

Rt  it,  lit ) = &  £  h  a“ 


r-i 


(19) 


We  are  now  interested  in  the  relationship  between  the  expected  value 
of  the  sample  covariance,  #■(&  Xij),  and  the  covariance 

*,-(*„  i 

a  , 

If  A,^  Ax,  B r~  because  the  measurement  noises 

are  spatially  uncorrelated. 


If  2,  -  A/ can  be  expanded  in  terms  of  the  multiplicative 
noise  and  the  image  corrupted  by  additive  and  Poisson  noise. 


The  mean  of  the  first  term  adds  to  the  variance  the  contribution  of 
the  additive  and  Poisson  noise.  It  is  equal  to 


E £(&/&•  Jj') -(fa))}-  +~  cU,).  (2i) 


The  mean  of  the  second  term  is  zero  because  /%?//„  J/  has  unity  mean 
and  is  statistically  independent  of  '  The  mean  of  the 

last  term  is  equal  to  the  product  of  the  2nd  order  non-central  moment 
°f  LdJ and  the  variance  of 


Therefore  we  obtain: 


)J  -  c(lt) 


(23) 


+  +  «*<) 


and  we  can  sum  up  the  cases  and  in  a  unique  formula 

The  difference  between  and  M  ^^Jcan  be  considered  as  a 
bias.  Once  this  bias  has  been  subtracted,  the  sample  covariance  (18) 
can  be  solved  for  the  object  using  the  two  broadband  dimensions  of 


For  the  assumed  case  of  isoplanatic  or  spatially  invariant  imaging, 
the  convolution  nature  of  (18)  allows  us  to  use  the  Fourier  transform 
to  rewrite  (18)  as: 

Sc'fe, jZx)  ~  0(Z,}@{ lx)  f—0-3.) .  (25) 

And  in  the  frequency  domain,  (24)  becomes: 
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EpJi.Ji)}  -  S;/ZlJli)*(i.i£)[<r.'Slil-QtI'I/iri  l26l 

<70  ° 

+  rpl.  -e)I\ P-/3 )  *  Si  (l  -e}  Z- *  jJJ/5 

If  the  noise  model  is  accurate  and  the  parameters  are  known,  the 
sample  covariance  can  be  corrected  for  the  effects  of  noise.  If  not, 
within  the  constraints  of  this  or  another  noise  model,  the  parameters 
can  be^estimated  and  used  to  correct  the  sample  covariance.  Let  us 
call  Zj)  the  corrected  sample  covariance. 

The  algorithm  previously  developed,  [1]  is  an  iterative  process  that 
uses  both  the  first  and  second  order  statistics  through  equations  (17) 
and  (25).  If  0  (2.)  denotes  the  X-th  estimate,  then  for  a 
particular  value  of  z(,  equations  corresponding  (17)  and  (25)  are 
written  as: 

6  (2,)  $(%)  =I(Z,)  1271 

&%,) Ss(i„zj]=  S-(z  jJ 

^  '28> 

=  <r2(ucr^)  SlZ.-ZJ 

0*0 

-  ff[I(Jr/l)Ik^)  *  SilZ,-£,  z, 

where  only  values  of  z^near  z_^are  used,  Jzf-Zj/<r0  •  Because  the 

sample  statistics  are  used,  these  equations  have  errors.  The 

equations  for  /z  -z  j>rc  are  dominated  by  errors  and  are  discarded. 

The  corresponding  components  of  are  not  computed,  greatly 

reducing  the  order  of  the  computation  (proportional  to  M 'rather  than 
y 

M  for  M  x  M  images).  These  equations  at  a  particular  value  of  z  ,  are 
-A/ 

solved  for  u  ( z )  in  a  least  squares  fashion  to  give 


2”  (ji)  ij ? i ) JA'  )S(zli  (z9) 

$  IZrZk  I  ^  rc 

Q  A)  ~  ~  /  7/  r  7i  77  ~ 

5Ll°  (zJ&ztJzj  +«ls(zti 

iz,-zj*r> 


A 


where  <*!  is  a  weighting  parameter  to  control  the  relative  importance  of 
the  sample  mean  and  sample  covariance  ( c<  =0  for  lz_J>  re  ).  These 
equations  are  solved  for  each  ^starting  at  the  origin  proceeding 
toward  the  diffraction  limit  in  a  spiraling  manner. 


The  parameters  of  the  noise  model  are  estimated  (if  desired)  after 
(29)  has  been  applied  to  all  ,z(.  The  estimate  is  obtained  by  a  least 
squares  solution  of 


+C,  Sd-lJ  i  cx  8(i,-zx)tC3l(zrzj 

for  the  coefficients  c0 ,  cx  ,  c2,  and  c^where  is  a 

precalculated  approximation  to  the  integral  of  (28).  The  coefficients 
c, ,  Cj  ,  and  c,  are  used  as  (/-><?£) ,  and  (1+^)1*/^  .  The 

square  root  of  ceis  used  to  scale  0l£,)  before  the  next  iteration. 

This  iteration  using  (29)  and  (30)  usually  converges  after  two  or 
three  passes.  Then  the  inverse  FFT  algorithm  is  used  to  compute  a 
reconstructed  object  image.  Finally,  by  varying  assumptions 
concerning  the  noise  processes,  we  obtain  several  reconstructions  and 
use  them  to  remove  the  noise  artifacts,  which  vary  between 
reconstructions. 


Finally,  our  signal-processing  method  for  restoring  atmospherically 
degraded  images  can  be  briefly  described  as  follows.  We  obtain  first 
and  second  order  statistics  of  the  degraded  images  from  40  to  300 
short-exposure  (about  10  ms)  images.  These  measured  statistics  are 
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then  related  to  the  underlying  object  and  atmospheric  characteristics 
by  a  series  of  integral  equations,  which  are  solved  by  a  combination 
of  "optimal"  signal -processing  techniques.  With  Fourier  methods,  we 
transform  the  integral  equations  into  a  set  of  non-linear  algebraic 
equations.  An  estimate  of  the  undegraded  object  is  then  obtained  by 
solving  linearized  subsets  of  these  equations  using  least  squares 
techniques.  The  restoration  technique  converges  quickly,  and  the 
required  computational  effort  is  roughly  proportional  to  the  number  of 
resolution  cells  in  the  image. 

Figure  1  summarizes  the  restoration  algorithm. 


» 
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Figure  1  -  Algorithm  for  reconstructing  the  object  image  from  sample 
statistics.  Once  the  iterative  process  converges,  a  reconstructed 
object  image  is  computed  with  the  inverse  Fourier  transform.  The 
final  reconstructed  image  is  obtained  after  several  preliminary 
reconstructions  enable  us  to  remove  the  noise  artifacts. 
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Section  2  -  Two  Approaches  for  a  New  Algorithm 
2-1  -  A  Recursi ve  Least-Squares  A1 gori thm 

Our  first  concept  for  a  recursive  algorithm  was  to  simply  transform 
the  previously  studied  least-squares  algorithm  into  a  recursive 
least-squares  algorithm. 


Let:  Ijzjbe  the  recorded  Nth  image  in  the  frequency  domain, 

©^Zjbe  the  ^th  estimate  of  the  object  obtained  with  N  images, 
77^/ and  (lLJ2)  be  the  first  and  second  order  sample 

statistics  computed  with  the  N  first  images. 


t 

A  test  of  convergence  is  conducted  on  the  sequence  Jind  for 

the  last  iteration,  Jmax,  the  estimate  )  is  obtained  (1  max 

A/  ■— 4 

may  depend  on  Z>  ).  When  the  next  image  (N+l)  is  available,  the 

estimate  of  the  object  is  initialized  for  each  z_.with, 

ctiij  =  (z  ).  Then  using  recursive  formulas  for  the  sample 

statistics,  the  equation  (29)  of  Section  1  can  be  rewritten  as 

follows: 


IZ,~  2, 1  £ 


The  first  drawback  of  this  real-time  algorithm  is  that  it  assumes  the 
object  and  the  statistical  characteristics  of  the  noise  processes  are 
stationary. 


Besides,  let  us  examine  recursive  formulas  for  the  sample  statistics 
of  the  recorded  images.  The  formulas  are  given  by: 


X  )  ~ a/*j_  J-  (’Zi.)  n * i  J-hii(Zc') 


(32) 
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A'k.i.i- ¥-A 


A 

A/t/  f  *</ 


V  «.,?,)  =  ,34) 

This  algorithm  would  require  the  storage,  at  each  step  N,  of  the 
quantities: 

I.  ]  for  each  z  . 

X M(g. )  for  each  /z^r 

0^/^  2  j  for  each  £.,z_.such  that  /z^ 

For  MxM  images ,  the  storage  requirements  would  be  about  10  to  20  times 
2 

M  and  would  make  the  use  of  a  small  computer  difficult. 

If  the  sample  statistics  did  not  have  to  be  stored,  a  small  computer 
could  be  used  more  easily.  This  would  require  storing  only  the 
current  estimate  of  the  object.  All  the  desired  information  in  the 
sample  statistics  is  ideally  contained  in  the  current  estimate  of  the 
object.  The  sample  statistics  recursions  (32)  and  (33)  would  be 
replaced  in  this  recursive  algorithm  by 
a  r 

- /£z  S<z<)  fzj  (35) 

This  allows  the  recursion  (31)  to  be  written  as  -J' 

el'h. i  -  ofuj f[fy<% Cl>  1  ^ n‘! * *rs<~‘ f]  <*> 

[jr,I  [efo  >  >3A  , /,  I  S'*.  >] 

A\ 

-?fTi h;  -tiffa  h  ■in*?,  a)  -  I  <*.  A  //  -■)  Jj 


l 


This  algorithm  is  similar  in  form  to  a  Kalman  Filter  where  Kalman 
weighting  functions  have  been  greatly  simplified.  A  Kalman  Filter 
approach  will  be  taken  to  obtain  a  more  general  result  with  a  better 
analytic  basis  for  its  form. 

2-2  -  A  New  Approach:  A  Kalman  Fi 1  ter  A1 gorithm 

The  problem  of  the  real-time  estimation  of  the  object  0(ZL),  in  the 
case  of  evolution  of  the  object  and  non-stationarity  of  the 
atmospheric  characteristics,  can  be  formulated  using  a  Kalman  Filter 
approach. 

As  for  the  existing  least-squares  algorithm,  the  assumption  of 
isoplanatic  imagery  will  enable  us  to  transform  integral  equations 
into  algebraic  equations  and  therefore  the  model  in  the  frequency 
domain  will  be  preferred  to  the  model  in  the  spatial  domain. 

The  remainder  of  this  paper  treats  the  theoretical  formulation  of  the 
Kalman  Filter  and  its  particularities  due  to  our  special  case. 


However,  before  studying  the  Kalman  algorithm,  a  slightly  modified 
notation,  taking  into  account  the  evolution  of  the  object  and  of  the 
statistical  characteristics  of  the  turbulence,  needs  to  be  introduced 

Let  the  Nth  image  measured  at  time  JN  be  represented  in  the 
spatial  frequency  domain  by 


Tn( Z(.)is  the  Fourier  transform  of  the  measured  image, 
is  the  optical  transfer  function, 

£>„(Z( ) i  s  the  Fourier  transform  of  the  object,  and 
n  ■  .  is  the  Fourier  transform  of  the  measurement  noise. 
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This  statement  is  illustrated  by  the  following  block-diagram. 


Spatial  Domain 


Frequency  Domain 
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the  arrow  symbolizes  the  two-dimensional  Fourier  transform 
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Section  3  -  Theoretical  Formulation  of  the  Kalman  Filter 

3-1  -  Definition  of  the  State  Vector  and  of  the  Observations  Vectors 

Both  the  real  part  and  imaginary  part  of  the  object  @(?t!  will  compose 
the  real  state  vector  X  «  which  we  want  to  estimate.  The  state 
vector  can  depend  on  time  (we  want  an  algorithm  that  permits  the 
evolution  of  the  object)  and  2k/  will  represent  the  state  at  time  N. 

At  each  time  N,  we  receive  a  new  image.  Let  us  call  the  Nth 

image,  in  the  frequency  domain,  and  at  the  frequency  •  This  image 
is  received  at  time  N.  With  this  image  we  can  define  a  vector  of 
observations  Y^.  Yu  is  composed  of  measurements  when 

fZ'-z^Hrc  ;  or  to  be  more  precise,  the  real  vector  )£.  is  composed  of 
the  real  and  imaginary  parts  of  ■!„  (z^) JT  fZA. 

3-2  -  Equations  of  the  Kalman  Fil ter 

The  purpose  of  the  Kalman  Filter  is  to  build  a  sequence  of  estimates 
of  the  states  N=l,...)  from  an  initial  estimate  1,0  of  the  state 

K,  •  Each  estimate  is  a  realization  of  a  random  process  and  is  then 
characterized  by  its  statistical  properties,  in  particular  the  mean 
and  covariance  matrix. 

Let  us  cal  1  X.,  ,/v.-<-»the  sequence  of  estimates,  resulting  from  our  Kalman 
Filter.  We  shall  require  for  our  estimate  X. ,  to  be  an  unbiased 
estimate  with  a  minimum  covariance  matrix 

In  order  to  build  X/y/^%  we  use  the  information  given  by  the  N  first 
images,  Xu/A/  is  then  called  the  a-posteriori  estimate.  The 
information  given  by  the  Nth  image  is  stored  in  the  vector  Y^  .  Using 
the  relation 

J*  ^c)  =  (38) 
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we  can  express  the  observations  Y#  as  a  quadratic  function  of  the 
state  X  : 

—A/ 

(3 

where  ^  is  a  deterministic  quadratic  function  and  l)A,  is  a  noise 
vector  of  zero  mean  and  covariance  matrix  R^.  This  is  the 
observation  equation  of  the  Kalman  Filter. 


The  observations  do  not  depend  linearly  on  the  state,  but  follow  the 
quadratic  equation  (39).  In  order  to  obtain  the  state  covariance 

we  have  to  use  a  Taylor  series  of  Ji(')  and  then  introduce  the 
Jacobian  matrix/^*  To  obtain  Q^as  a  function  of  only  » 

KV) and  we  have  to  restrict  the  Taylor  series  to  the  first  order 
items. 


As  shown  in  the  following  Section  3-3,  the  noise  of  observations  is 
not  Gaussian.  However,  we  will  use  only  the  first  and  second-order 
moments. 

These  "non-linear  observations"  and  the  use  of  the  Jacobian  matrix  to 
linearize  them  about  the  current  estimate  results  in  an  algorithm 
known  as  the  extended  Kalman  Filter.  This  solution  of  will  be  optimum 
in  the  sense  of  minimum  covariance,  but  not  in  the  sense  of  maximum 
likelihood. 

Estimation  of  the  state  given  the  first  image,  X// 


If  the  a-priori  state  estimate,  7C^  ,  is  unbiased,  and  if  we  look  for 
an  unbiased  estimate  of  X//  as  a  linear  combination  of  {f  and  7^  »  we 
can  prove  that  7  must  follow  the  equation: 


2f,.  =  Ju.  *K,(X-  %(lu 
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The  gainK,,  also  called  Kalman  gain,  will  be  chosen  such  that  the 
covariance  matrix  of  ,  is  minimum. 

Let  us  call  Ci/A  the  covariance  matrix  of  2fy0  •  Then,  expanding  (40) 


r‘/o 

about  yC.  gives 

— Vo 


This  relation  permits  us  to  calculate  Ch  as  a  function  of  Ci/A , 


(41) 


K/,  and  Kt  . 


J 


Cl/t  -  [J-KH, )  Q (I- KH)  +  K •  K, K  (42) 

Qy  is  a  definite  positive  symmetric  matrix  associated  with  a  definite 
positive  quadratic  form  that  will  be  minimum  for 

»i«>  D,- H,  C,aHt-iR 

and  we  shall  have: 


n 

i 


(43) 


C»  =  Q  - «,  H,  4 

Prediction  of  the  state  at  the  second  image  ^  from  X 


(44) 


'// 

The  evolution  of  the  state  between  the  times  N  and  N+l  is  given  by  an 
equation  of  evolution  of  the  following  type: 


2^vy/  -  J-  (ZN  jUn+i)  *  K+ 


(45) 


where  ^is  a  deterministic  function  that  can  house  a  very  general 
form,  where  U^is  an  input  vector  and  where  i^,?/is  a  noise  vector  of 
zero  mean  and  covariance  matrix  QNfl> 


» 
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To  simplify  the  presentation,  we  suppose -^/linear: 

y  4,, 


(46) 


Using  the  estimate  which  is  the  "best  estimate"  in  the  sense  of 
minimum  covariance  available  for  _/T  ,  we  can  predict  the  unbiased 
es  ti  ma  te  2^  for  2^  by : 


Tv, 


c 


z/i 


is  called  the  a-priori  estimate  of_^,  its  covariance  matrix, 
,  is  given  by 


G,  -  4  C( 


y. 


h,  f  Q, 


(48) 


Using  the  second  image  and  repeating  the  operation  we  performed  to 
obtain  from  ^  we  obtain  a-posteriori  estimate  of  :  2y  and 
i  ts  covari  ance  C?/z  • 


Using  alternate  filtering  and  evolution  equations,  we  can  compute 
(2f */*;  N=1 ,...).  The  Kalman  Filter  equations  can  be  summarized  by 
the  following  set  of  equations: 
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Evolution  Equati ons 


'a  21'  ^ 


(50) 


3-3  -  Particularities  of  the  Noise  of  Observations  ^  y)  and  Form  of 
the  Equation  of  Observation 

Another  particul ari ty  of  this  Kalman  Filter  is  the  model  of  the  noise 
of  observations  ft*.  This  noise  results  from  a  combination  of  noises 
of  different  sources: 

Noise  due  to  the  turbulence  of  the  atmosphere  which  is  modeled  by 
a  random  transfer  function  of  the  telescope,  5^  at  the  frequency  Z- 

Noise  of  measurement 
noise). 

In  order  to  understand  the  noise  process,  we  can  use  a  block-diagram 
and  some  simplified  notation. 

Given  a  time  N,  we  can  omit  the  index  N  appearing  in  the  notation  and 
we  define  K;  and yj  as  follows: 


(51) 

(52) 


h  (additive,  multiplicative  and  Poisson 

't./v 
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J3T  is  the  image  received  by  the  telescope  at  the  frequency  J/.  The 
real  and  imaginary  parts  of  y  are  elements  of  the  state  vector, 
y t  is  the  measured  image  at  the  frequency  Z,- . 

By  definition  of  the  transfer  function  5^  at  the  frequency  zt-  we  have 

y'  =  5}  Xt-  (53) 


and  by  definition  of  the  noise  of  measurement 

/ 


y;  -  X-  +  *1;  i 

Lji  =  +  rl;  < 


(54) 


(55) 


To  help  understand  the  form  of  the  observation  equation  and  the 
particularities  of  the  noise  ''V,  it  is  interesting  to  calculate 


where  E^pjdenotes  the  expected  value  of  the  random  process  p. 


The  real  and  imaginary  parts  of  the  terms  E 


will  compose  the 
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elements  of  the  deterministic  vector  Also,  the  real  and 

imaginary  parts  of 

-X 

noise  vector  Y\  . 

Therefore,  by  construction  will  be  deterministic  and  the  noise 

vector  77^  will  have  a  zero  mean. 


{viruiE1 


)  will  compose  the  elements  of  the 


Calculation  of  i  form  of  the  function^  and  matrix 


The  two  random  processes,  noise  of  turbulence  and  noise  of 
measurements,  are  statistically  independent.  Also,  the  noise  of 
measurement  has  zero  mean.  This  allows  E To  be  simplified, 

r}*)  (68) 

The  parameters  of  the  complex  function  are  composed  of  the 

second-order  moments  of  the  noises  of  turbulence  and  measurement: 

Eft-pand  t(wp. 

Knowing  that  E  /X  is  a  real  number,  the  elements  of  will 

have  one  of  the  following  forms: 


IU k  nf)tb, s, V •L(f(yl,r\  V- (toUtit >)) ■!-(%) ■J-q>pv }7' 4/% q Vj 


or 


Then,  the  derivatives  of  these  expressions  with  respect  to  the  real  or 
imaginary  parts  of  ?c  or will  give  the  elements  of  the  Jacobian 
matrix/^.  These  elements  will  be  one  of  the  following  types: 


R. ctylffi.S/Jj z#e  <V  £ (S‘  <*■ 


0. 


The  second-order  moments  of  the  noise  of  turbulence  are  the  parameters 
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of  the  Jacobian  matrix//^. 


The  Noise  h  and  its  Covariance  Matrix  R  , 

- A - '  yi/ 

-V. 

The  elements  of  are  given  by  the  real  and  imaginary  parts  of 

ir  ttjfisf-ia,#* 

Some  statistical  properties  of  the  noise  n.,  will  be  given  for  the 
covariance  matrix  of  which  is  involved  in  the  Kalman 

equations. 

The  elements  of  RK.  will  be  of  one  of  the  three  following  forms: 

L  far.  if -tWfl) &/%£-£ j))  = 

j  ?/?*)] 

E(l~  Y y*  -ty •  %>*) ) h*  *  4  n 

-  k  a\  -(\t.  )}\%  y/i  -  ;  -t/ :4)l 

The  general  formula  for  Cov(y,y;  )RR  ^  is  deve1oPecl  in  Appendix  1 
and  will  permit  us  to  calculate  the  matrix  , 

The  general  formula  for  Cov(y,  y,*  f  is 

Crv(p  =  *. 

fit)  -  sRm/if) 

+  */*  L(s!sf  I  FfWI  -  1,  *  Ffs.  U)  £ /rf/l  ) 


^  7  f- ^  "ft1* )  *  *4  Ff'it  )F  fa,  t\  V  > 

+£(wji\/  C )  -  =  ■  %  nj 1  /  (y  a;  ) 


U  - 


I 
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The  following  table  gives  the  orders  of  various  moments  in  the  terms 
of  Covt^^l. 


For  the  noise  and  for  the  noise 

of  turbulence  of  measurement 


in  terms  of  the  1st  line  4 

in  terms  of  2nd  and  3rd  lines  2 

in  terms  of  the  4th  and  5th  lines  1 

in  terms  of  6th  line  no  moment 


no  moment 
2 

3 

4 


We  can  notice  first  that  the  sum  of  the  orders  of  the  moments  of  the 
noise  and  turbulence  and  of  the  noise  of  measurement  is  always  equal 
to  four,  and  second  that  third-order  moments  of  the  noise  of 
turbulence  have  disappeared  because  they  were  multiplied  by  the 
first-order  moment  of  the  noise  of  measurement,  that  is  by  zero. 

Conclusion 


The  calculation  of  the  second-order  moments  of  the  noise  of 
observation  nK  will  require  the  knowledge  of  the  1st,  2nd,  3rd  and 
4th  order  moments  of  the  noises  of  turbulence  and  measurement.  This 
is  due  to  the  nonlinearity  of  the  observation  equations. 

3-4  -  Particularities  of  the  Covariance  Matrix  of  the  State  Vector 

One  particulari ty  of  the  state  vector  comes  from  the  fact  that  it 
contains  real  and  imaginary  parts  of  complex  numbers;  and  this 
particularity  will  confer  a  special  structure  to  its  covariance 
matrix. 

Let  us  consider  "T  the  state  vector  at  the  time  N  and  an  estimate  of 
it  7^..  The  covariance  is  given  by 
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a.  v*.  -zj 


(61) 


4  is  equal  to  Q  if  \sh/K  or  t0  Iv-i/t  ,'  iis  comP°sed  of 

the  real  and  imaginary  parts  of  U  ■)  and  we  can  also  define  the 
complex  vector  ^whose  elements  are  <4  U,).  An  estimate  of^is  given 
by^  .  The  dimension  of  ^  and  is  twice  the  dimension  of  £  and 


We  can  also  define  two  complex  matrices  and  ,^as  follows: 

Q,  =  E\(&- iV  XA  -&)Tj  ( «  > 

c4  -  E  If i  y£  -%)*] 


If  p  and  q  are  some  elements  of^,,  then  £>  will  be  composed  of  terms 
of  the  form  Cov(p,q)  and  will  be  composed  of  terms  of  the  form 
Cov(p,q*~) . 

Also,  is  composed  of  terms  of  the  following  form: 

Cn-f/LtDjLfi))  i Re (63) 
C  Cl-  fr  (f)j  2c  f ) )  -  i  2c  ft -riff,  9  )  _  2l -(p ,  f  JJ 

C«r  )^-J 


We  want  to  express  Ch  as  a  function  of  £T  and  &  . 


Then,  assume  that 


(64) 
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(C/j/is  equal  to  the  block-matrix, 


Using  the  sets  of  formulas  (62)  and  (63)  we  can  show  that, 

A,  -  x  e,jc„  ■  £„  ] 

Because  X^^^is  a  symmetric  matrix  and  is  an  antisymmetric 

matrix,  we  can  rewrite  (68)  and  (66)  as 


(66) 

(67) 

(68) 


_  T  P  r  "V'  — 

B,  = 


and 


/">  _  /  .A  /'C-C 

Wt'  "  I - 


*\S 


Separating  the  terms  concerning  ^  and  we  obtain 

AJ  M 


(69) 


(70) 
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Section  4  -  Particular  problems  of  the  application  of  a  Kalman  Filter 
to  Speckle  Imaging 

4-1  -  Decomposition  of  our  Large-Scale  System  into  Subsystems 

When  implementing  a  Kalman  algorithm  to  treat  speckle  images,  the 
first  problem  we  meet  is  a  problem  of  dimension. 

A  typical  dimension  of  a  digitized  image  has  been  256  by  256  elements 
(pixels).  We  will  use  the  discrete  Fourier  transforms  of  the  measured 
images  as  input  data  at  a  discrete  set  of  two-dimensional  frequencies 
Z‘X).  Using  the  symmetry  properties  of  the  Fourier  transform 
and  the  aperture  limits  will  reduce  the  set  of  frequencies  that  must 
be  treated.  The  result  of  estimating  the  object,  in  the  frequency 
domain,  is  also  a  discrete  set  of  values  of  the  Fourier  transform  of 
the  object.  This  set  will  have  the  same  dimension  as  the  one 
concerning  the  image  and  will  be  treated  by  the  same  FFT  algorithm  to 
restore  the  estimate  of  the  object  in  the  spatial  domain.  Then,  the 
state  of  our  Kalman  Filter,  ,  will  typically  have  about  51.5DD 
elements  (256  x  256  %7fj 4). 

This  is  a  large-scale  system  that  will  require  a  specific  algorithm  of 
decomposition.  However,  the  decomposition  of  the  whole  system  into 
several  subsystems  of  reasonable  size  will  come  naturally  when 
examining  the  structure  of  the  observations.  As  a  matter  of  fact,  the 
measurement  yt-yy*  does  not  bring  any  useful  information  whenever 

.  Thus,  when  we  want  to  estimate  a  particular  */,  we  shall 
select  a  subsystem,  from  the  full  system  composed  of  all  the 

/  if.  t 

7'S  such  that  y^yy  is  one  of  the  measurements  used  in  the  full 
system. 

Then  the  state  vector  of  the  subsystem  J?  is  a  vector  restricted  to 

,  *>■ 

the  real  and  imaginary  parts  of  the  T^Ssuch  that  /z.-zy^r„  and  the 
observations  vector  is  composed  of  all  the  measurements  y,'y*  that  can 
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be  formed  with  the  X’a  and  Xj*  belonging  to  the  state  vector  of  the 
subsystem. 

Since  £  is  very  small  compared  to  the  aperture  R,  the  dimension  of  the 
state  of  the  subsystem  is  very  small  compared  to  the  dimension  of  the 
state  of  the  full  system. 


This  decomposition  is  actually  an  approximation.  Let  us  go  back  to 
the  full  system  and  let  us  consider  three  frequencies  Z- »  2y  and 
such  that:  fZi  '2jMtc,  !Zj-2j/<%,  and 

Then,  the  two  measurements  y,  ui  and  will  belong  to  the  vector  of 
observations  of  the  full  system,  but  not  y ^  iH. 

j/L. 

However,  the  measurement  will  be  helpful  for  the  estimation  of 

'V  and  this  estimation  of  Z,'  will  affect  the  estimation  of  2V,  by  the 
intermediary  of  the  observation  Therefore,  ,  which  is  an 

observation  of  the  subsystem ,  will  affect  the  estimation  of 
7 l  even  if  /z.-Z/is  greater  than  r_.  Thus,  the  subsystems^/’  and 
Jf'  are  not  independent  even  if  /Z -Zyyfr 

*-k 

This  shows  that  the  Kalman  filter  propogates  the  information  given  by 
a  particular  measurement  The  whole  set  of  frequencies 

Z>  and  that  the  estimation  of  will  be  affected  by  the  whole  set  of 
measurements 

whatever  Z>  ,  7  ,  and  2/  may  be.  Obviously,  the  dependence  of  the 

r  ™  4^ 

estimation  of  7;  upon  the  measurements  (  )JA/~//2"'  will  be 

more  or  less  important  depending  on  the  distances  /Z£  ~Z  -/and 
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The  Effect  of  Propagation  of  the  Information  on  the  Covariance  Matrix 
of  the  State  of  the  Kalman  Filter 


If  we  initialize  fusing  the  first  image  and  the  information  this 
image  gives  with  the  measurements  ^-^-such  that  }2t  r^,  then 

the  covariance  matrix  of  the  state  will  be  initialized  by  a  matrix 
with  a  band  structure  and  the  width  of  the  band  will  depend  on  the 
number  of  frequencies  inside  the  circle  of  radius  ra. 

As  a  result  of  the  phenomena  of  propagation,  if  we  run  the  Kalman 
Filter,  the  width  of  the  band  of  the  covariance  matrix  will  increase 
with  each  new  image  and  some  covariance  terms  will  appear  between  two 
estimates  T  and  such  that 

Implication  on  the  Choice  of  the  Subsystems 

The  algorithm  on  the  subsystems  must  restore  this  propagation  of 
the  information  in  order  to  use  each  observation  as  much  as  possible. 
To  achieve  this  goal,  we  must  make  an  adequate  choice  for  the 
subsystems. 

Assume  that  we  make  the  following  choice  of  subsystems.  We  partition 
the  set  of  frequencies  _z^into  "circles"  of  radius  re  such  that  each 
respectively  Xt  )  belongs  to  one  and  only  one  "circle"  (respectively 
subsystem).  The  treatment  of  one  subsystem  associated  with  a 

subset  of  frequencies,  will  provide  the  estimates  of  all  the  7j'  s  such 
that  .  Then  the  treatment  of  all  the  subsystems  will 

provide  an  estimate  of  the  whole  object. 


Nevertheless,  with  this  choice  of  subsystems,  we  would  not  restore  the 
phenomena  of  propagation  of  the  information  along  the  set  of 
frequencies  because  we  would  make  two  subsystems  appear  independent, 
whereas  they  are  not  in  the  full  system.  Also,  the  use  of 
measurements  would  not  be  optimal:  the  estimation  of  ,7.  belonging  to 


t 


* 


r 


» 


i 


t 


i 


* 


t 


» 


belonging  to^f?!c  ,  such  that  ^  ■  is  different  from 
zj/^  ,  would  not  use  the  observation  tjj  ^  .  This 
discarding  of  the  measurements  that  correspond  to  crossproducts  of 
state  elements  in  separate  subsystems  is  undesirable.  Because  such  a 
choice  is  inappropiate,  we  must  chose  a  sequence  of  subsystems  that 
retains  these  measurements.  There  are  several  closely  related 
viewpoints  which  allow  the  use  of  the  needed  measurements. 

First,  we  can  choose  a  sequence  of  K  subsystems </(  1 ) ,</( 2 ) , . . .  ( K ) 

such  that  two  adjacent  subsystems^m)  and^m+l)  have  a  common  part 
in  their  state  vector.  Then  the  part  of  the  object,  which  is  common 
to^/im)  and.J'Wl)  and  estimated  by  the  Kalman  filter  on  the  subsystem 
yf(m),  will  be  used  by  the  Kalman  filter  on^/(m+l).  Thus,  using  the 
a-posteriori  estimate  fo r  we  shall  improve  the  a-priori  estimate 

for  J/m'lJ 


J  and  offy 

but  kr. 


Second,  we  can  expand  the  observations  vector  used  with  a  given 
subsystem  J..  to  include  all  measurements  that  involve  at  least  one 
element  of  the  state  vector.  Then,  all  of  the  data  will  be  used. 

Finally,  we  can  view  all  of  these  subsystem  constructions  as 
particular  cases  of  treating  the  observations  of  the  full  system  in 
subgroups  and  restricting  the  elements  of  the  state  vector  being 
updated.  The  updates  are  then  treated  in  a  sequence  which  treats  all 
of  the  observations  and  all  of  the  state  vector.  The  subgroups  of 
either  the  observations  or  the  state  vector  or  both  must  overlap. 

This  overlapping  characteristic  of  the  algorithm  will  propagate  the 
information  on  the  estimation  of  the  state.  Thus,  the  state  estimates 
will  be  interdependent  and  correlated.  However,  the  corresponding 
terms  of  the  covariance  will  still  be  zero  and  the  norm  of  the 
covariance  will  be  larger  than  computed  by  the  full  system. 

let  us  illustrate  the  first  viewpoint  by  a  one-dimensional  example. 


I 
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Let  us  consider  a  set  of  Fourier  transforms  of  one-dimensional  images. 
The  set  of  frequencies  has  been  discretized  into  M  discrete 
frequencies  separated  by  intervals  of  constant  length  ;z0=Q,  z,,... 
im.r  We  choose  the  subsystems  so  that  the  state  of  the  subystem 
^(2+1)  is  obtained  from  the  state  of  the  subsystem^/d)  by  dropping 
the  object  at  the  frequency  z1_l ,  and  inserting  the  object  at  the 
frequency 

The  treatment  of</*(l)  will  provide  an  estimate  of 

The  treatment  of  will  provide  an  estimate  of 

The  treatment  will  provide  an  estimate  of  ~XX)  -/' 

After  the  treatment  of</(l),  we  keep  the  estimate  of  x„because  it  has 

been  made  using  all  the  information  jfe  ‘jo^i  .•••  We  a1so 

use  the  estimates  of  x(,  x^,...  x„.,  as  the  a-priori  estimates  for  the 

Kalman  Filter  of ^(2).  The  estimate  of  x.,  we  obtained  withal)  did 

*•  ' 

not  use  the  information  given  by  y,  . 

Then,  with  J?(2)  we  shall  again  estimate  x, ,  now  using  The 

result  will  contain  information  about  y<-^/*given  through  the  a-priori 
estimation  of  x, ,  and  information  about  , ...  given 

through  the  observations  vector  of<P(2).  We  can  repeat  this  operation 
until  the  Kalman  filter  has  been  applied  to  the  last  subsystem.  We 
start  again  with  the  treatment  of <^R  1 )  when  we  receive  a  new  image. 
Such  an  algorithm  will  use  the  important  measurement  information  and 
will  propagate  the  information  along  the  set  of  frequencies. 

If  the  discretization  of  the  Fourier  transforms  of  the  images  is  very 
fine,  corresponding  to  a  large  field  of  view,  we  shall  have  too  many 
points  inside  a  circle  of  radius  r0  and  the  state  of  any  subsystem  may 
still  be  too  large  to  be  treated  easily  and  quickly  by  a 
mini -computer.  In  this  case,  we  can  either  reduce  the  field  of  view 
and  treat  subsections  of  available  images,  or  formulate  the  Kalman 
filter  using  a  more  restricted  update  of  the  state  vector. 


-35- 


4-2  -  The  Kalman  Gain  Matrix 

Computing  the  Kalman  gain  matrix  usually  requires  most  of  the 
computational  effort  associated  with  the  filtering  equations.  A 
straight  forward  treatment  of  the  Kalman  gain  equation  for  the  full 
system  would  require  the  inversion  of  an  enormous  matrix  whose  order 
is  the  same  as  the  observations  vector.  Fortunately,  by  using  the 
decompostion  into  subsystems,  the  order  of  the  matrix  that  needs  to  be 
inverted  is  greatly  reduced.  However,  because  of  the  high  number  of 
images  and  subsystems,  the  computation  effort  required  for  inversion 
would  still  be  quite  large. 

The  required  computational  effort  can  be  reduced  further  if  we  forego 
some  of  the  generality  of  the  Kalman  formulation.  We  shall  assume 
that  our  system  has  no  dynamics.  By  "no  dynamics",  we  mean  that  the 
equation  of  evolution  ^  (X ,  UAt/  )4>  becomes  7^  -J*  •  This 

simplification  restricts  the  class  of  objects  being  estimated  to 
constant  images  described  by  the  general  diffraction  limited  Fourier 
transform.  Although  we  do  not  allow  any  evolution  for  the  object,  we 
allow  the  statistical  characteristics  of  the  noises  of  turbulence  and 
measurement  to  change.  For  example,  changes  in  the  seeing  conditions 
and  measurement  conditions  are  quite  coitwon. 

How  much  the  required  computational  effort  can  be  reduced  depends  on 
the  degree  to  which  the  Kalman  gain  calculation  can  be  simplified. 

Two  alternatives  will  be  explored  briefly.  First,  the  Kalman  gain 
calculation  will  be  approximated  by  a  difference  equation.  Second, 
further  simplifying  assumptions  will  be  used  to  allow  the  Kalman  gain 
to  be  precomputed  to  various  degrees. 

Simplification  of  the  Notation 

As  the  system  has  no  dynamics,  the  evolution  equations  of  the  Kalman 
Filter  become, 
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for  the  state: 


I 


t 


► 


Vv  l72) 

for  its  covariance: 

AS -t  /  Jfj  d-^/sS  (73) 

ihe  a-priori  variables  at  time  N+l  are  equal  to  the  a-posteriori 
variables  at  time  N.  Therefore,  we  can  get  rid  of  the  double  index 
notation  and  j^and  (^would  now  represent  the  a-posteriori  variables, 
which  result  from  the  Kalman  filter  at  time  N.  The  word  estimate,  if 
nothing  else  is  precisely  expressed,  will  refer  to  the  a-posteriori 
estimate. 


With  the  new  notation,  the  filtering  equations  become: 
Gain: 


State : 


Covariance: 


=  c,.,  hJd;' 

D/I/  -  (bU iCv-i^a/ 

L  J  2 


) 


2,  -  2»->  *  K*-  (  y»  - 

1  c  --  c«-,  c*-> 


(74) 


(75) 

(76) 


Development  of  a  Difference  Equation  for  the  Kalman  Gain 

Let  us  now  consider  the  equations  of  the  gain  giving  and  : 


(77) 

178) 
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And,  in  order  to  develop  a  difference  equation  for  the  gain, 
define  the  following  matrices: 

let  us 

f 

Gt/  - 

(79) 

G+/  ~  Fsv+I 

m 

! 

d/j  (-•'As—  t  (((a/ 

dd  (D„,  )  ~  Da/*i 

(81) 

(82) 

and  (respectively  and  are  computed  with  the  a-priori 
estimates  of  the  state  times  N+l  and  N,  2^  and Ttien> 
^/^(respectively  A  )  will  depend  on  the  variation  of  the  state 
estimate  between  times  N-l  and  N,  but  it  can  also  depend  on  the 
variation  of  the  noise  moments  between  times  N  and  N+l  as  these 
moments  are  the  parameters  of  the  Jacobian  matrix  (respecti vely  the 
covariance  matrix/?). 

We  must  notice  the  minus  sign  in  (81).  The  Kalman  Filter  causes  a 
decrease  in  the  covariance  matrix  of  the  state  estimate  at  each 
iteration. 

As  we  shall  see  later,  the  last  one  among  these  variation  matrices, 

can  be  approximated  in  the  first  order  by  a  linear  function  of 
the  matrices  >  AR*,,  and  AC^. 

The  equation  (78)  can  now  be  rewritten  as  follows: 

Kw  =  (&-,  -A  G  )( G/  4v/ ))  ( 83 ) 

and  assuming  small  variations,  if  we  only  Keep  the  first  order  terms, 
we  obtain 

» 
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K,„  -  L  -(& 

This  is  the  first  step  to  the  difference  equation  we  are  looking  for, 
and  we  now  need  to  develop  the  expression  ofA(/P4/T/)  as  a  function  of 
AHm  tAtfm,  and AQ,-  Expanding  VKfj  and  keeping  only  the  first  order 
terms,  we  obtain: 

Dfii\  ~ Ov(l -  Ov  ( C/l'Hv  ~A^,Ca  Hu 


If  we  denote  by  AN.  the  following  matrix: 

A/ 


-  4 LC*  Hl c*-' ^  ^  ^ y/  ^ 


we  obtain 


If  we  use  the  following  approximation: 


-  A  // 
^  '  U. 


7  -  a  ft 

-1-  ^  dv  /  ^ 


we  obtain 


MC, )  -  4 


The  difference  equation  for  the  Kalman  gain  becomes 


/C  -  ■* /-  ^  *  6-,  (90) 


ol[H»  0^-4C/4>>^9U 

-) 

This  recursive  formula  will  permit  calculating  K,„  from  /%  , 
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p 


anc*Zj  f\Art/. 


W 


-/  -/ 

Then  assuming  we  know  Z?v  •  we  will  not  have  to  calculate  DAr/  by  a 

matrix  inversion  in  order  to  obtain  A^y/.  With  the  Kalman  gain  A^we 
can  easily  get  the  state  estimate  ,  its  covariance  the 

variation  matrices,  and  finally  0. '  • 

v’y/ 


Precomputed  Kalman  Gain  Matrices 

Finally,  we  will  examine  conditions  under  which  the  Kalman  Gain  matrix 
and  state  covariance  can  be  partially  or  fully  precomputed.  The 
Kalman  gain  and  state  covariance  depend  on  each  other,  but  only  the 
Kalman  gain  is  required  for  the  state  estimation  recursion.  First,  we 
will  determine  the  conditions  which  reduce  the  Kalman  filter  to  the 
recursive  least  squares  algorithm  of  Section  2.  Then,  we  will 
consider  the  effect  of  removing  some  of  these  simplifying  conditions. 

The  recursive  least  squares  algorithm  (36)  updates  the  state  at  only 
one  frequency  and  uses  all  the  observations  within  rQ  of  that 
frequency.  Under  these  conditions,  the  Kalman  filter  has  a  two 
element  state  vector  and  a  larger  observation  vector.  The  Kalman  gain 
and  state  covariance  equations  can  be  rewritten  in  the  form  where  the 
matrix  inversion  lemma  has  not  been  applied, 

K  =[C,  ’HYhJ'h'k'  (92) 

C,j  - <  ///?  ’/■/] 


(93) 


Assuming  that  the  atmospheric  and  measurment  statistics  are  stationary 
gives  constants  for  the  covariance  R  and  the  moments  used  in 
calculating  observation  function?^)  and  its  Jacobian  Assuming, 
that  the  Taylor  series  from  which  the  Jacobian  is  obtained  is  expanded 
about  the  correct  result  rather  than  the  current  estimate  yields  a 
constant  Hy.  Under  these  conditions,  the  state  covariance  and  Kalman 
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gain  equations  can  be  solved  to  yield 

C„-  7 

Kv  -  ^Ih1r,h]'h1R' 


The  resulting  recursion  becomes 

r^(0A  /it)]  -A'JI  K 

!_7*  ( eA  (lL )  J  /■£’*>  (&*-•(- 1  h 


(96) 


where  the  current  estimate  can  be  used  to  evaluate  the  Kalman  gain 
expression. 


This  result  is  directly  analogous  to  the  least  squares  recursion  (36) 
of  Section  2,  but  further  simplification  is  required  to  complete  the 
comparison.  The  difference  of  observations  and  their  predictions  are 
the  same  for  both  recursions, 

4  5/7, ,  2 j  )=  X-  3  ■/  w)  ( 97 1 


Each  row  of  the  Jacobian  matrix  is  the  partial  derivative  of  one  of 
the  observations  with  respect  to  the  state  and  has  one  of  three  forms: 

-i 


u  for  ft  > 

2.  for 

3.  for  Imty/]; 

If  the  observations  are 
variance,  R  is  given  by 


[&[Ov5sa^)}7„  (iC/u  /?„?/ 
Jf*  %  <2,1/  Z*  K  <?,  f  17,7/ 1 

assumed  to  be  independent  with  the  same 


(98) 


Now  the  Kalman  gain  matrix  of  (96)  is  quite  close  to  the  gain  matrix 
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of  (36).  The  recursion  derived  from  (97)  and  written  in  a  form 
similar  to  (36)  is 

j)l  JFUteji;))]  +  r  uT^sc 

-  _J  L- 


W-  fS!]' 

^'WlRe(eU))W iii)&ui),B-fai*zJSsaJZ)f 


A=  1  (Z;,z,-)  , 

1 

S=I  /eC^) S^z,yz j) . 

Izrijl4‘r' 

Note  that  (36)  used  a  diagonal  matrix  inversely  proportional  to  B  for 
the  covariance  of  the  state  C^which  provides  the  same  normalization 
for  both  the  real  and  imaginary  parts  of  the  estimate  and  does  not 
cross  couple  them.  Also,  we  have  not  treated  lower  frequencies, 
J_z.j<r,  ,  in  (99). 

Increasing  the  dimension  of  the  state  vector  does  not  substantially 
change  the  recursion  (96).  The  Kalman  gain  can  still  be  derived 
a-priori  and  evaluated  using  the  current  estimate.  If  we  assume  that 
the  observation  statistics  vary  only  by  a  multiplicative  constant,  we 
only  have  to  compute  a  scalar  weighting  factor  to  replace  the  1/N 
factor  in  (96). 

If  we  allow  the  observation  statistics  to  vary  more  generally,  the 
Kalman  gain  cannot  be  derived  a-priori.  However,  if  we  assume  some 


1 


restrictive  description  of  how  the  statistics  vary,  we  can  approximate 
the  Kalman  gain  a-priori.  In  particular,  assume  that  the  observation 
statistics  are  slowly  varying.  For  example,  the  seeing  conditions  as 
described  by  the  mean  optical  transfer  function  SU)  or  rp  ,  could  be 
slowly  varying.  Then,  SU)  could  be  estimated  over  intervals  less 
than  the  decorrelation  time  of  the  seeing  conditions  and  used  in  the 
recursion  (95)  that  assumes  stationari ty. 

Conclusion 

We  have  developed  a  Kalman  formulation  for  the  Speckle  Imaging 
problem.  The  intent  was  to  develop  a  formulation  that  would  allow  the 
use  of  mini -computers  for  Speckle  Imaging.  This  was  accomplished  by 
using  decompositions  of  the  full  system  into  manageable  subsystems. 
These  subsystems  result  from  natural  decompositions  based  on  the  very 
limited  area  of  coherence  in  the  aperture  caused  by  atmospheric 
turbulence.  Thus,  the  limitation  of  nature  which  requires  that 
Speckle  Imaging  is  needed  aids  in  making  the  numeric  solution 
manageable. 


Although  this  formulation  was  developed  under  very  general  conditions 
for  the  image  and  disturbance  models,  we  have  examined  simplifications 
and  approximations  which  reduce  the  computational  requirements.  These 
simplified  Kalman  filters  approach  a  recursive  version  of  the  least 
squares  algorithm  previously  developed. 


Appendix  1  -  Calculation  of  the  Terms  used  in  the  Formation  of  /\4/, 
the  Covariance  Matrix  of  the  Noise  of  Observations 

A  -  Calculation  of  the  variance  and  covariance  terms  of  the  real  and 
imaginary  parts  of  two  complex  random  processes 

If  Ur)  is  the  expected  value  of  a  real  random  process  r,  we  naturally 
define  the  expected  value  of  a  complex  random  process  by: 

E(p)  =  U  Ref p)  )  +  j  E(  Imfp)  ) 

Then,  we  can  define  the  complex  covariance  of  two  complex  random 
processes  p  and  q  by: 

Cov(p,q)  =  E(  (p  -  E ( p)  )  (  (q  -  E(q)  )  ) 

Using  to  the  definition  of  E(-),  we  can  commute  the  operation  E(-) 
with  any  of  the  following  ones:  Re(*),  Im(-),  and  conjugation. 
Therefore,  we  can  develop  the  two  following  calculations: 

1  >  Gm,-  (Prf)  =  Qnr(. AVfV tj  ReUCf ) )  ~  (in  (p)y[ m (cj) ) 

Cnr (f?4 £  (p)j  L,  STL'  (iWRjfe  (</)) 

2)  C<nr  (fj  f) 

(<j) )+jC<nr(j*  (p)j 

Combining  Cnp^cjj  and  <?)  we  ^ind  the  three  following 

formulas: 

CnJ ~  J  Rt^Cnrfjij  cj)  -t  Cci-fa  <fjj 
Gtrl-J Ke(p ~  J  Xm  JjTn-Cfttf)*  Cov-(/°j  R j] 
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We  shall  apply  these  formulas  to  p  =  y^,  y,  and  q  =  y^  and  develop 
below  a  general  formula  for: 


and  then 


Cov(p.q)  =  Cov(  y.y*  ,  y  y*  ) 
1  i  £  l 


Cov(p.q')  =  Cov(  y- y*  ,  y^y*  ) 


To  obtain  Cov(p,q)  from  Cov(p.q),  we  shall  interchange  the  two  indices 
l  and  X  • 

B  "  “Complex  Covariance"  of  the  Two  Random  Complex  Processes  y  y -  and 


We  want  to  calculate 


Covl  y.y*  ,  fyyf  )  «  E(y.  -  E(yf.  yp  E(y^  y' 


*  x, 


*m  +  ^  + 


Expanding  pq  =  y;  y.*y,  y*  gives 
A  X 


for  m  =  i ,  j ,  k,  or  1 . 


Pi  =  A- *; f  ^  ^ 
+  *4s<$(Wtf+njxAsl) 

+x£ss;(i\:*X-<-kS) 


rH  1  N‘77 


X  ,'X  K 


+ i  ty  x<^  )(r\n  x^5  *  ^  ^  )  • 


To  calculate  E(pq),  we  take  the  sum  of  the  expected  value  of  each  term 
of  the  previous  sum  and  we  use  the  statistically  independent  of  the 
noise  of  turbulence  and  the  noise  of  measurement.  Also,  because  of 
the  zero-mean  of  the  noise  of  measurement,  the  two  underlined  terms 
have  a  zero  mean.  Then,  the  moments  of  p  and  q  are: 


Combining  these  equations  we  obtain  the  general  "complex  covariance": 
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ABSTRACT 

*  ... 

This  report  deals  with  astronomical  imaging  through  a  turbu¬ 
lent  atmosphere  when  the  commonly-made  assumption  of  isoplanicity 
is  not  valid.  Integral  equations  are  developed  which  relate  the 
object  to  the  mean  and  autocorrelation  function  of  a  series  of 
short-exposure  images.  Useful  approximations  to  certain  functions 
describing  the  effects  of  the  atmosphere  are  developed  and  tested 
using  extensive  numerical  computations.  These  approximations  lead 
to  a  considerable  simplification  in  the  analysis  of  the  effects  of 
nonisoplanatic  imaging. 


I 
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1.  Introduction 


The  purpose  of  this  report  is  to  describe  the  effects  of 
atmospheric  turbulence  on  astronomical  imaging  when  the  object  is 
so  large  that  the  imaging  is  not  isoplanatic.  We  begin  by  review¬ 
ing  the  equations  which  describe  optical  imaging  under  the  usual 
conditions  of  Fraunhofer  diffraction.  We  then  develop  some  equa¬ 
tions  which  relate  the  object  to  the  first  and  second  moments  of 
the  image.  Under  certain  conditions  which  are  frequently  encoun¬ 
tered  in  practice,  these  equations  are  well  conditioned  and  can  be 
solved  to  yield  a  reasonable  estimate  of  the  object.  The  imaging 
equations  have  been  developed  by  Sherman  [1]  and  this  report  is  an 
extension  of  his  analysis. 

Unless  further  assumptions  and  approximations  are  made,  the 
imaging  equations  will  involve  complicated  functions  and  time- 
consuming  numerical  integration.  A  key  contribution  of  the 
research  described  in  this  report  is  the  development  of  useful  and 
accurate  approximations  which  allow  us  to  obtain  closed-form  ana¬ 
lytical  results  which  describe  the  effects  of  nonisoplanat ic  imag¬ 
ing.  Results  obtained  using  the  approximate  model  are  compared 
with  the  corresponding  results  obtained  using  numerical  integration 
in  order  to  assess  the  accuracy  of  the  approximation.  As  we  shall 
demonstrate,  the  approximate  model  yields  very  accurate  results 
over  a  wide  and  useful  range  of  parameter  values.  Furthermore, 
bounds  on  the  spatial  frequencies  for  which  the  model  is  valid  are 
established,  and  these  bounds  are  related  to  physical  phenomena. 

The  imaging  configuration  with  which  we  are  concerned  is 


shown  in  Fig.  1.  When  the  object  is  incoherent,  the  image  inten¬ 
sity  distribution  is  given  by  the  superposition  integral 

CO 

l(x)  =  /  /  0(y)  s(x,  y)  d^y  (1) 

—CO 

where 

x  is  the  angular  position  (in  radians)  in  the  image 

plane, 

y  is  the  angular  position  (in  radians)  in  the  object 

plane, 

l(x)  is  the  image  intensity  distribution, 

0(y)  is  the  object  intensity  distribution, 

and 

s(x,y)  is  the  point-spread  function  of  the  atmosphere- 
telescope  system. 

We  assume  that  the  conditions  for  Fraunhofer  diffraction 
apply  [2].  Under  these  conditions,  the  complex  amplitudes  in  the 
object  and  aperture  planes,  and  those  in  the  aperture  and  image 
planes,  are  related  by  Fourier  transforms.  The  Green's  function 
for  a  point  source  at  angular  position  y  in  the  object  plane  is 
given  by 

00 

U(x,  y)  =  A  /  /  F(z,  y)  exp  {jkz  •  (y  -  x)}  d^z  (2) 

—00 

where 

z  is  the  aperture  position  in  meters, 
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A  is  a  scale  factor, 

U(x,y)  is  the  complex  amplitude  at  angular  position  x  in 
the  image  plane  resulting  from  a  point  source  at  y 
in  the  object  plane, 
k  is  the  wave  number, 


F(z,  y)  is  the  complex  amplitude  at  z  in  the  aperture  plane 
resulting  from  a  point  source  at  y  in  the  object 
plane. 

We  write  F(z,  y)  in  the  form 


F(z,  y)  =  \|)(z,  y)  a(z)  exp  (j9(z)} 


where 


<|>(z,  y)  is  the  complex  transmittance  function  of  the  atmo¬ 
sphere, 

a(z)  is  the  pupil  function  of  the  aperture, 


0(z)  is  the  phase  aberration  function  of  the  lens. 

The  transmittance  function  <Kz,  y)  can  be  expressed  in  the  form 


t(i(z,  y)  *  exp  {x(zi  y)  +  j4>(z,  y )} 


where 


x(z,  y)  is  the  log  amplitude  disturbance 


and 


$(z,  y)  is  the  phase  disturbance  due  to  the  atmosphere. 

In  terms  of  the  quantities  defined  above,  the  incoherent 
point-spread  function  becomes  [2] 

s ( x,  y)  =  U(x,  y)  U*(x,  y) 

CD  00 

=  I Al 2  /  /  /  /  F(B  +  w,  y)  F*( 8 ,  y)  (5) 

—00  —00 

2  2 

•  exp  {-jkw  ♦  (x  -  y)}  d  8  d  u>. 

F,  and  therefore  s,  are  of  course  sample  functions  of  random  pro¬ 
cesses  which  must  be  described  in  terms  of  their  statistical  prop¬ 
erties.  We  turn  in  the  next  section  to  the  development  of  expres¬ 
sions  which  relate  the  first  and  second  moments  of  the  image  inten¬ 
sity  distribution  to  the  object  intensity  distribution  and  to  the 
statistical  properties  of  the  atmosphere. 

2 .  First-  and  Second-Order  Image  Statistics 

Sherman  [1]  has  derived  expressions  for  the  mean  and  the 
autocorrelation  function  of  the  image  intensity  distribution  under 
the  assumption  that  x(*>  •)  and  $(•,  ♦)  are  Gaussian  random  pro¬ 
cesses.  Similar  results  have  been  obtained  by  Fante  [3]  and 
others.  We  briefly  review  these  results  in  this  section.  Detailed 
derivations  are  contained  in  Appendix  A. 

Taking  expected  values  of  both  sides  of  (1)  yields 

00 

E{l(x) }  =  J  /  0(y)  E{s(x,  y)}  d2y.  (6) 
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The  mean  point-spread  function  has  been  shown  to  be  [1,  3] 


00 

E(s(x,  y)}  =  1 A | 2  /  /  exp  {-D(w,  0)/2}  •  exp  {-jku  •  (x  -  y)} 

—CO 


•  /  /  a (0  +  w)  a(8)  exp  (j6(3  +  »)  -  j 6 ( B )}  d2w  d20. 

—  00 

(7) 

where  D(A,  6)  is  the  two-point  wave-structure  function  whose  prop¬ 
erties  will  be  discussed  in  some  detail  in  section  3. 

The  autocorrelation  function  of  the  image  intensity  distribu¬ 
tion  is  given  by 


E{l(x1)l(x2)}  =  //  //  0(yi)0(y2)  •  E{s(x1>yi)  s(x2>y2)}  &2 y  ^ y  ^ 

— OO  —00 

(8) 


where 


E{s(x1,y1)s(x2,y2)}  =  |a|A  //  //  //  //  ,0^  ,  ,u>2 ) 

—00  —00  —00  —00 

•  M(si,“1.y1»62,w2,y2)  exp  {-jk<D1  •  (X]l  -  y^  -  jku>2(x2 

,2  .2.  .2  .2 

•  d  6j  d  82  d  tij  d  »2< 


The  quantities  P(  )  and  M(  )  in  (9)  are  defined  by 


?{Zl,ul,&z,(a2)  =  a(61  +  ^  )  a^)  a(B2)  a[82  -  o>2) 

•  exp  {je(e1  +  +  i©( b2 )  -  je(e2  -  u>2 )}  do) 


and 
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>  y^>  ^2 *  t°2 ’  y2^ 


exp  <-T 


j  Dfwj.  0)  +  D(a)2,  O) 
j  +  d[61  -  32  +  ^  +  w2,  h(yi  -  y2)] 

2  +  d[Bx  -  02,h(yi  -  y2)]  -  D[B1  -  B2  +  «»1,h(y1 

D[0!  "  B2  +  W2’  h(yl  “  y2^ 


The  derivation  of  (9)  to  (11)  neglects  terms  involving  the 
autocorrelation  of  x(*»  *)  and  the  cross-correlat ion  of  x(*»  *)  and 
4> C *  ,  •).  Fante  [3]  has  argued  that  the  impact  of  these  terms  is 
small.  For  further  details  on  these  derivations,  see  Appendix  A. 

p(B1 ,  u)^ ,  B2>  uj2)  describes  the  effects  of  the  finite  aper¬ 
ture  and  of  lens  aberrations  on  the  second-order  statistics  of  the 
image,  while  m((^,  w^,  y^ ,  B^ ,  w2 ’  descrihes  the  effects  of 

atmospheric  turbulence.  Sherman  [2]  has  argued  that  lens  aberra¬ 
tions  can  be  neglected  if  turbulence  is  strongest  near  the  ground 
in  the  sense  that 


Ah  h'+Ah 

f  CZ(h)  dh  >  /  C  (h)  dh 

0  °  h’  n 


where  C  (h)  is  the  refractive  index  structure  constant  at  height 
n 

2 

h.  We  will  assume  condition  (12)  to  be  valid  for  all  profiles 
of  interest,  and  will  therefore  neglect  the  effects  of  aberrations 
in  the  remainder  of  the  report. 
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We  introduce  the  following  changes  of  variables: 


8+6  6  -  B 

1  2  1  2 
&  =  ;  A  6  = 


-  yl  +  y2  .  yl  '  y2 
y - - -  ;  Ay  =  - ? 


to  -  U)  (0,  +  to 

-12.  12 

to  =  - - -  ;  Ato  =  - - - 

2  ’  2 


Making  these  variable  changes  and  neglecting  lens  aberrations  in 
(10)  and  (11)  yield  the  equations 


P'(8,  AB,  to,  Ato)  =  a(B  +  AB  +  w  +  Ato)  a(B  +  A3)  a(B  -  AB) 


•  a(B  -  AB  +  w  -  Ato) 


M'  (to ,  Ato,  AB,  Ay)  = 


/  D(to  +  Ato,  0)  +  D(to  -  Ato,  0)  N 

exp  <  -  j  +  D(2AB  +  2Ato,  2hAy)  +  D(2AB,  2hAy) 

\-  D(2AB  +  Ato  +  to,  2hAy)  -  D(2AB  +  Ato  -  to,  2hAy) 


We  now  turn  to  a  more  detailed  discussion  of  the  two-point 
wave-structure  function  D(*,  •). 
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3.  The  Two-Point  Source  Wave-Structure  Function 

Kon  and  Feizulin  [4]  have  derived  the  following  expression 
for  the  two-point  source  wave-structure  function: 


D(A,  6) 


=  0.033 


r|i),2kV5/3  "  c2(h> 

6  /  m  q  n 


dh 


(16) 


where 

A  is  the  distance  between  test  points  in  the  aperture 

plane  (meters), 

6  is  the  distance  between  source  points  in  the  object 

plane  (meters), 

k  is  the  electromagnetic  wave  number, 

Km  is  5.91/1q,  where  1Q  is  the  inner-scale  size  of  the 

turbulence , 

H  is  the  height  of  the  turbulent  medium, 

2 

C  (h)  is  the  refractive  index  structure  constant  at  height  h, 
n 

L  is  the  distance  from  object  to  aperture  plane, 

is  Rummer's  function  (a  confluent  hypergeometric 
f unct ion) . 

Rummer's  function  is  the  confluent  hyper geometrc  function  [3] 


jFj (a,  b,  z) 


az  (a)(a  +  1)  z^  (a)(a  +  1  )(a  +  2)  z^ 

b  (b)(b  +  1)  2!  +  (b )  (b  +T)(b  +  2)  3!  + 

(17) 
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r 


i 


> 


i 
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This  is  a  very  messy  function  to  work  with,  and  we  seek  reasonable 
approximations.  In  the  literature  [6],  it  has  been  common  practice 
to  develop  approximate  expressions  for  ^ F j  C a ,  b,  z)  for  the 
limiting  cases  |z|  <<  1  and  |zl  >>  1.  We  observe  that  for  | z 1  << 

1,  we  may  neglect  the  higher  order  terms: 

jF^a,  b,  z)  -  1  +  ,  |  z  1  «  1,  (18) 

jFjU,  b,  z)  -  1  -  ~  ,  I z  1  «  1.  (19) 

For  |z|  »  1,  we  employ  the  asymptotic  expression  [5]  valid  for 

Re [ z ]  <  0: 

jF^a,  b,  z)  n  rTb"-^)'  (~z)~a»  (z|  »  l-  (20) 

Since,  in  the  problem  at  hand,  a  <  0  and  z  <  0,  jFj(a,  b,  z)  is  an 
increasing  function  of  z.  For  |z(  sufficiently  large  we  have 

jF^a,  b,  z)  -  1  "  j^-;--)  (~z)~3  ,  |z|  »  1.  (21) 

Fante  [3]  has  used  these  limiting  values  of  Rummer’s  function  to 
derive  approximations  to  (16).  His  approximate  expressions  are 

(H* 


A  + 


H 


dh 


(22) 


O  1  *  o 

D(A,  6)  -  k 2  j  CZ(h)  d 
i  n 


where 


I 
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I 


d(p)  = 


3.25  lo'1/3  IpI2, 


2.92  |p 


5/3 


IpI  «  0.17  l 

O 

IpI  »  0.17  l 

O 


(23) 


Upon  introducing  the  imaging  problem  in  section  1,  we  used 
angular  rather  than  linear  units.  To  reflect  this  in  the  structure 
function,  we  define  the  angular  position  vector  6'  as 


6' 


(24) 


We  also  assume  that  the  objects  are  astronomical  and  hence  that  L 
>>  H.  With  this  assumption  and  the  change  of  variables  shown 
above,  we  find  that 

9  H  9 

D( A,  hi’)  «  k  /  C  (h)  d( A  +  h6 ’ )  dh.  (25) 

0  n 

A  key  issue  is  the  behavior  of  d(p)  when  p  ■  0.17  1  .  Thus 

o 

far  we  have  only  developed  approximations  valid  in  the  regions 

|p|  «  0.17  1  and  |p|  »  0.17  1  .  It  is  not  known  at  this  point 
o  o 

which  of  the  two  approximations  models  the  behavior  of  d(p)  in  the 
region  of  uncertainty.  Indeed,  it  is  not  known  if  either  approxi¬ 
mation  is  valid  in  this  intermediate  region.  To  determine  more 
closely  the  behavior  of  D(A,  h$'),  b»  z )  was  computed 

numerically  by  summing  terms  in  the  hypergeometr ic  series.  The 
result  was  compared  to  the  quadratic  approximation  and  to  the  5/3- 
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power  approximation  ((19)  and  (21),  respectively).  The  results  of 
the  numerical  analysis  of  Rummer's  function  are  detailed  in  Appen¬ 
dix  B.  In  brief,  it  was  found  that  the  quadratic  approximation 
(19)  begins  to  diverge  rapidly  beyond  the  breakpoint  p  =  0.17  1  , 
while  the  5/3-power  approximation  closely  follows  Rummer's  function 
over  the  entire  range  tested.  We  conclude  that  the  quadratic  law 
is  valid  only  for  very  small  arguments  and  that  the  5/3-power 
approximation  is  reasonable  for  all  z.  Henceforth,  we  will  use  the 
5/3-power  approximation  of  d(p)  for  all  values  of  p. 

With  an  analytical  expression  for  the  wave  structure  func¬ 
tion,  we  are  ready  to  consider  the  concept  of  the  atmospheric  cut¬ 
off  frequency  .  We  associate  with  the  spatial  frequencies 
that  can  be  recovered  from  the  blurred  image  by  long-term  averag¬ 
ing.  The  cutoff  is  determined  by  the  first-order  statistics  and 
is  generally  defined  [7]  as  the  value  of  u)  that  results  in  an  expo¬ 
nent  of  -1  in  the  function  exp  [-1/2  D(u>,  0)].  That  is, 

exp  [-1  d(u>c,  0)]  =  exp  [-1].  (26) 


Solving  for  algebraically,  we  obtain 


u 

c 


2 

2.92  k2  vi 

o 


3/5 


(27) 


where  we  have  defined  y  to  be  the  zero-order  moment  of  the  refrac- 

o 

tive  index  structure  constant 
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H  2 

M  =  /  C  (h)  dh. 

°  o  n 


In  general,  we  define  the  -order  moment  y  to  be 

n 


y  E  /  hn  Cvh)  dh. 
n  0 


These  moments  will  be  useful  in  later  developments. 


4.  Imaging  Equations 

Substituting  the  expression  for  the  structure  function  given 
by  (25)  into  (15),  we  obtain  the  following  expression  for  the 
mutual  coherence  function  M'(w,  Aw,  AS,  Ay): 


M’(ai,  Aaj,  AS,  Ay)  = 


(d(u+Ao))  +  d(u-Au>) 

+  d(2Ag+2hAy+2Au>)  +  d(2AS+2hAy) 

_ 

-  d(  2AS+2hAy+A(u-uj)  -  d  ( 2AS+2hAy+Aa>-Ki)) 


dh  . 


We  make  the  following  observations  about  the  mutual  coherence  func¬ 


tion: 


-  2  2  2 

1.  For  frequencies  satisfying  |w)  +  |Au>|  <  ,  the  mutual 

coherence  is  essentially  independent  of  AS  and  Ay,  and 
the  entire  aperture  appears  to  be  spatially  coherent. 


That  is, 


« 


r 


M'(w,  Aw,  A3,  Ay)  »  Aw). 


2.  For  spatial  frequencies  satisfying  |Aw|  >>  , 

|w|  »  w^,  the  mutual  coherence  is  nearly  zero: 


M'(w,  Aw,  A3,  Ay)  *  0. 


3.  For  spatial  frequencies  satisfying  |Aw|  < 

|w|  »  there  is  nonzero  coherence  in  the  region 


| 2A3  +  Aw|  <  u  and 

c 

1  Ay  I  < 

Ay  .  In 

c 

this 

case 

the 

aperture  is  spatially 

coherent  in  regions 

the 

size 

of  W^,  provided  that 

the  two 

source 

points 

are 

not 

separated  by  a  distance 

greater 

than  Ay  . 

The  ability  to  recover  the  high  spatial  frequency  information 

from  the  second-order  statistics  of  a  set  of  speckle  images  depends 

upon  the  nature  of  the  mutual  coherence  function  in  case  3  (i.e., 

|w|  >>  u  ,  |Ao)|  <  a)  ).  In  the  remainder  of  this  report  we  will 
c  c 

concentrate  on  analyzing  the  coherence  function  under  these  condi¬ 
tions. 

We  will  now  consider  the  effects  of  the  turbulence  on  the 
image.  As  discussed  previously,  we  neglect  lens  aberrations.  We 
have  shown  that 

OB  OB 

E{l(x1)l(x2}}  =  //  //  0(y1)0(y2)E{s(x1,y1)s(x2,y2)}  d  yjd  y*  (8) 

"  —OB  —OB 

where 
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I 


* 


oo  co  oo  oo 


eMv  yx)six2>  y2)}  =  4>Al4  //  //  //  //  p,(w»  Aui»  AS> 


—  OO  —00  —00  —  oo 


_  2—  2 
•  M'(u),  &o),  AH,  Ay)  d  &  d  A3  exp  [-jkw^ 


2  2 

(X1  “  yI^  "  Jki)2  *  ^X2  "  y2^d  "l  d  V 


(9) 


P’(w,  Aw,  8,  A8)  =  a(8  +  A0  +  w  —  Aw)  a(3  +  A3)  a(8  -  A3) 
•  a(3  -  A3  +  w  -  Aw) , 


(14) 


and 


M'(w,  Aw,  A3,  Ay)  = 


exp 


d(w+Aw)  +  d(w-Aw) 

2  H 

~  /  C2(h)  j  +  d(2AB+2hAy+2Aw)  +  d(2A8+2hAy)  | dh 

0  °  . 

y-  d(2A8+2hAy+Aw-w)  -  d(2A3+2hAy+Aw+w)> 


.  (30) 


To  simplify  the  notation  we  introduce  the  function 


N(w,  Aw,  Ay)  =  //  //  M'(w,  Aw,  A3,  Ay)  P’ (w,  Aw,  3,  A3)  d20  d2A3 

—  OO  — oo 

(31) 


The  integration  over  A3  effectively  samples  the  function 
P'(w,  Aw,  8,  AB)  at  A3  =  -Aw/2.  Hence 


} 


♦ 
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Y 


N(u,  Au,  Ay)  *>  //  P'(u),  Au,  B,  AB  =  - 


Tf  <’* 


//  M'  (u,  Au,  AB,  Ay)  d2AB. 


(32) 


For  convenience  we  define 


F4(w,  Au)  h  JJ  P'^u,  Au,  $,  AB  =  -  d2B 


=  //a(B  +  u+  -^jafB-^|jafB  +  2 


Au 


a  |  B  +  u  -  ~  )  d2B 


(33) 


and 


Mq(u,  Aw,  Ay)  =  //  M'(u,  Au,  AB,  Ay)  dZAB.  (34) 


From  (8),  (9),  (14),  and  (31)-(34),  we  find  that 


E{l(x1)l(x2)}  =  4|A|4  //  //  0(yi)  0(y2)  //  //  Mo(u,  Au,  Ay) 


—00  —00 


•  F4(w,  Au)  expj-jk^  •  (xj  -  yL)  -  jku2 

*  (x2  “  *2^  d2“l  d2“2  d2yl  d2y2’  <35> 


This  equation  is  valid  for  |u|  »  u  ,  and  lAu|  <  u 


We  may  express  0(y)  and  Mq(w,  Au,  Au)  as  inverse  Fourier 
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transforms : 


0(y)  = 


00 

//  O(w')  exp 
—  00 


[- jfcw ' 


y]  d2w’, 


(36) 


M0(“»  Aa)’  W  =  Ui 


//  M  (to ,  Aw,  Aw’)  exp  [-jkAw’  •  Ay]  d  Aw1. 

— OO 

(37) 


Substituting  these  expressions  for  0(y)  and  Mq(w,  Aw,  Ay)  into  (35) 
and  simplifying  yields 


lUC-iW-j))  ■  4IM4  ^  /*  //  //  5  U  -■5f)  0  U 

(k)  -®  -®  -»  \  '  \ 


Aw1 


“  -  2 
•  M  (u j,  Aw,  Aw’)  d  Aw*  F.  (w,  Aw) 
o’’  4  ’ 


exp  [~jkto1 


*l  ~  3^2 


2  2 
d  u)^  d 


(38) 


Or,  taking  the  Fourier  transform  of  both  sides  of  this  equation, 


Ffe[r(x1)rCx2)]}  h  w  (a>v  o>2) 


e  41a! 


4  (  2tQA 
(k)2 


•  M  (u>.  Aw,  Aw1)  F,  (w,  Aw)  d  Aw’ 
o’*  4 


(39) 


> 
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From  (38)  and  (39)  we  conclude  that  the  effect  of  nonisoplanat ism 
is  to  smooth  the  second-order  statistics  in  the  Aw'  plane.  The 
behavior  of  the  smoothing  function  Mq(w,  Aw,  Aw')  in  the  Aw'  plane 
determines  the  feasibility  of  solving  (38)  for  the  Fourier  trans¬ 
form  of  the  object  0(w).  If 

Aw'  <  w  ,  (40) 

c  c 

then  the  integral  equation  (38)  is  well  conditioned. 

In  the  remainder  of  this  report  we  shall  consider  the  impact 
of  the  mutual  coherence  function  M'(w,  Aw,  AB,  Ay)  and  the  smooth¬ 
ing  function  M  (w,  Aw  Aw1)  on  the  second-order  statistics.  Some 
o 

key  issues  include: 

1.  Is  it  possible  to  find  simple  analytical  expressions  for 
M’(w,  Aw,  AB,  Ay)  and  Mq(w,  Aw,  Aw’)? 

2 

2.  How  does  the  random  process  C  (h)  affect 

r  n 

M'(w,  Aw,  AB,  Ay)?  What  are  the  pertinent  properties 
2 

of  C^(h)  that  impact  the  imaging? 

3.  What  generalizations  can  be  made  about  the  impact  of  the 
smoothing  function  Mq(w,  Aw,  Aw')  on  the  imaging? 

5.  The  Effects  on  Nonisoplanat ism 

In  this  section  we  will  explore  the  effects  of  nonisoplana- 
tisra  by  deriving  reasonable  closed  form  approximations  for 
M  (w,  Aw,  Ay)  and  M  (w,  Aw,  Aw').  We  shall  temporarily  assume  that 


(41) 


|w|  »  0)^  ,  |Au>|  <  ,  1 2 A0  +  Acu |  <  Wc ,  |Ay!  <  y^ ; 

these  conditions  will  be  explored  in  more  detail  later.  We  begin 
by  reviewing  (30): 

M'U,  tin,  A0,  Ay)  = 


exp 


.2  H 

*2  I  Cn(h) 
0 


d^(w  +  Aw)  +  d^(w 


Aw) 


+  d^(2AB  +  2hAy  +  2Au>)  +  d^(2Af5  +  2hAy)  ]  dh  > 
d  ,  (  2Ag+2hAy+Aw-w)  -  d  ,  (  2A8+2hAy  +  Aw+u))y 


where 


(30) 


d(p)  =  2.92  |p 


5/3 


(23) 


We  have  numbered  the  d(»)  functions  with  subscripts  to  expedite  the 
discussion. 

The  conditions  set  forth  in  (41)  insure  that  |u)|  is  large 


while 

|Afl|, 

|Aw|,  and  |Ay| 

are 

small. 

For 

the 

moment  we  will 

assume 

that 

1  to |  dominates 

the 

other 

terms 

in 

the  arguments  of 

d^(*),  d^*),  d^(»),  and  d^(»).  This  assumption  is  suspect  since 
the  term  Ay  is  always  multiplied  by  h,  which  may  take  on  large 
values  (2  •  10  J  in  the  range  of  integration  [0,  H] .  However,  we 
shall  tacitly  accept  the  assumption  that  |  In  Ay  1  <  |w|  at  this  stage. 
Later  on,  we  shall  explore  the  limitations  and  consequences  of  this 
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assumpt ion. 


With  the  assumption  that  |to|  dominates  the  other  terms  in  the 
arguments  of  d^(*),  d^C*),  d^(*),  and  d^(*),  we  may  expand  each  of 
these  functions  in  a  power  series  by  employing  the  binomial  expan¬ 
sion  [5] 


(1  +  x)  =  1  +  a  x  + 


a(a 


2! 


1)  2  .  a(a  -  1) ( a  -  2)  3  . 

x  +  3j  x 


I  x|  <  1 . 


(42) 


For  the  problem  at  hand, 


d(to  +  6)  =  2.92  |to  -*■  5|5/3 

=  2.92  [(to  +{)•(«+  6)]5/6 


=  2.92 


,  ,  26  •  to  6*5 

to  |  1  +  — - —  +  - - — 

0)  •  0)  to  •  to 


3/6 


(43) 


As  before,  the  dot  denotes  inner  product.  We  have  assumed  for  each 
function  that  |to|  >  |S|,  so  we  may  expand  the  term  in  parenthesis 
in  a  binomial  series  with 


25  •  id  5  •  5 

x  =  - - +  — - — 

tO  *  to  to  •  <0 


The  resulting  expression  for  d(to  +  5)  is 


(44) 


d(5  +  5)  =  2.92  ( 


to 


- 15/6 
to) 


1  + 


55 


to  t  5 


3u 


(|  "  if  cos20) 


(45) 
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where 


0^  =  angle  between  w  and  Aw 

©2  =  angle  between  w  and  -Aw 

9^  =  angle  between  -w  and  (2A8  +  2hAy  +  Aw) 

0^  =  angle  between  +w  and  (2A0  +  2hAy  +  Aw). 

We  note  that 


91  = 


•n  —  0 , 


0 


3 


and  from  the  trigonometric  identity  cos  (tt  -  0)  =  -cos  0  ,  we 
conclude  that 


2  2 
cos  0^  =  cos  ©2 


2  2 
cos  0„  -  cos  0, 
3  4 


We  next  observe  that  the  mutual  coherence  function  M'  (w, 
Aw,  A8 ,  Ay)  will  be  maximized  for  fixed  w  and  Aw  when  (2A8  + 
2hAy)  "  -Aw.  This  leads  us  to  conclude  that  over  the  region  where 
M'  (w,  Aw,  A0,  Ay)  is  most  significant,  the  terra  (2A8  +  2hAy  +  Aw) 
will  be  nearly  colinear  with  Aw.  Hence: 


2 

cos  0. 


2 

cos  0, 


2  2 
cos  0.  =  cos  0. 
3  4 
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In  the  remainder  of  this  report  we  shall  employ  the  approximation 
above,  and  we  will  refer  to  8^  as  simply  9. 

Unfortunately,  we  cannot  expand  d^(»)  and  d^(»)  by  the  bino¬ 
mial  expansion  since  we  cannot  guarantee  that  any  term  in  the  argu¬ 
ment  will  dominate  for  either  function.  However,  we  recognize 
that  d^(*)  and  d^(*)  will  only  profoundly  impact  the  shape  of  the 
mutual  coherence  function  M'  when  the  magnitude  of  the  argument  is 
close  to  Hence,  we  will  approximate  d^(p)  and  d^(p)  by  a  sim¬ 

pler  function  that  is  a  good  fit  in  some  sense  when  p  B  w  .  The 

c 

choice  for  the  approximating  function  is 


d(p)  «  c.p  •  p  +  c 
1  o 


where  c^  and  cQ  are  constants  selected  to  "optimize"  the  fit  for 
|p I  “  Wc.  The  motivation  for  this  quadratic  approximating  function 
will  become  apparent  shortly.  Applying  this  approximation  to 
d^(*)  and  d^(*),  we  have 


d3(2A8  +  2hAy  +  2Aw)  «  c^AB  +  2hAy  +  2Au>)  •  (2AB  +  2hAy  +  2Aw) 


+  V 


d4(2AB  +  2hAy)  -  c^AB  +  2hAy)  •  ( 2AB  +  2hAy)  +  c  .  (53) 


To  simplify  the  notation  we  introduce  the  following  func- 

t ions : 
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f^w,  Aw,  A0,  Ay)  =  d^(w  +  Aw)  +  d2(w  -  Aw) 

-  d_(2A0  +  2hAy  +  Aw  -  w)  -  d, (2A0  +  2hAy  +  Aw  +  w) ,  (54) 

D  D 


f 2 ( Aw,  A0,  Ay)  =  d3(2A0  +■  2hAy  +  2Aw)  +  d^(2A0  +  2hAy).  (55) 


With  this  notation  we  can  express  (30)  as 


f  k2  2 

M'(w,  Aw,  A0,  Ay)  =  exp  < - r  /  C  (h) 

1  o  n 


*  [fjC". 


Aw,  A8,  Ay) 


+ 


(Aw,  AB,  Ay)]  dh 


) 

I 


(56) 

Combining  (47)-(50)  and  (52)— (53)  with  (54)  and  (55),  we  obtain 

f^(w,  Aw,  AB,  Ay)  "  ~19.47(l  -  cos^  0)  (w  •  w)  ^  [ A3  •  AB  + 

2 

(2hAy  +  Aw)  •  AB  +  hAw  •  Ay  +  h  Ay  •  Ay],  (57) 


f2(Aw,  AB,  Ay)  »  [8AB  *A3  +  (16hAy  +  8Aw)  •  AB  +  4Aw  •  Aw 

+  8hAw  •  Ay  +  8h^Ay  •  Ay]  +  2cq.  (58) 

Substituting  these  approximate  expressions  for  f^(w,  Aw,  AB,  Ay) 
and  f2 (Aw,  A0,  Ay)  into  (56),  we  obtain 
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I 


t 


I 


M'(w,  Aw,  AB,  Ay)  « 


exp 


/cj[8AB  •  AB  +  (16hAy  +  8Aw)  •  AB  \ 

k2  2 

f  2  1  l 

+  4Aw  •  Aw  +  8hAw  •  Ay  +  8h  Ay  •  AyJ  +  2cq 
i  -  19.47  1  -  -j  c os 20  (w  •  w)"1/6[A3  •  AB 

dh  ‘ 

\  +  (2hAy  +  Aw)  •  AB  +  hAw  *  Ay  +  h^Ay  •  Ay]y 

(59) 

Employing  the  definition  of  the  moments  of  the  refractive  index 

2 

structure  constant  C^(h)  given  in  (29),  we  may  perform  the  integra¬ 
tion  with  respect  to  h.  Computing  the  integral  and  collecting 
terms  involving  AB  •  AB  and  AB  yields 


M'(w,  Aw,  AB,  Ay)  *» 


exp 


8c.  y  -  19.47[l  -  4-  cos2  e)(u>  •  w)_1/6  y  (  AB  •  AB 
1  o  V  3  /  o 


16c^y^Ay+8c^PoAo>-19.47^1-  cos20^(w*w)  ^^(2PjAy 


+UQ Aw) 


•A0+  *c,  [8y„ AyAy+Uv  Aw*Aw+8y,  Aw*Ay1  t  2c  |) 

'  1 v  2  ^  o  1  oc 


-19.47^1-  -j  cos20^(w*w)  ^^(y^AwAy+y^Ay’Ay) 


(60) 


To  make  this  seemingly  complicated  expression  more  tractable, 
we  will  make  the  following  definitions: 


rt^  =  4CjPo  “  9.73  |^1  -  — •  cos2  ©j  (w  •  w)  1  ^  yQ,  (61) 


-  26  - 


(62) 


^2  =  8CjP^Ay  +  4c^UqAw  -  9.73  j  cos^  0)  (w  •  w) 

•  (2p^ Ay  +  PqAw), 

Tl~  =  c.  (4p  Ay  •  Ay  +  4p  Aw  •  Ay  +  2p  Aw  •  Aw)  +  c  p 
3  11  1  o  oo 

-  9.73  ^1  -  ~  cos^  6^  (w  •  w)  (p^Aw  •  Ay  +  p^Ay  •  Ay). 

(63) 

We  note  that  n2  is  a  vector  while  and  n^  are  scalars.  With 
these  definitions  (60)  can  be  expressed  very  simply  as 

r  j  ) 

Aw,  AB,  Ay)  »  exp  j -k  [AS  •  AB  +  n2  •  AB/^  +  ] j . 

(64) 

Since  this  expression  involves  only  quadratic  and  linear  terms 
in  AB,  we  may  perform  the  integration  with  respect  to  AB  called  for 
in  (34).  Hence,  with  the  approximations  we  have  made,  we  may  ob¬ 
tain  an  approximate  analytical  expression  for  Mq(w,  Aw,  Ay).  This 
was  the  prime  motive  for  approximating  the  d(»)  functions  by  qua¬ 
dratic  polynomials. 

Completing  the  square  in  the  exponent  of  (64)  and  performing 
the  integration  with  respect  to  AB,  we  derive  the  following  approx¬ 
imation  of  Mq(w,  Aw,  Ay): 
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Mo(id,  Aw,  Ay) 


CO 

J  I  M'(w,  Aw,  AS,  Ay)  d2A8 

—CO 


~T~~  exp  {-k2(n3  -  n2  •  njAnJ}.  (65) 

r  ni 


The  term  (n^  ~  tl2  •  r^/Ah^)  can  be  evaluated  by  direct  substitution 
of  (61 )— ( 63) .  The  result  is 


(n3  -  n2  •  hjAnJ  =  Yl  ( \  -  wJ/uo)  Ay  •  Ay  +  Y-^Aw  »  Aw  +  coiio, 

(66) 

where  we  have  made  the  following  definitions  to  maintain  the  sim¬ 
plicity  of  the  notation: 


Finally, 


M 

o 


'l  =~4Cl  -  9.73  (l 


1  2  .  -a-1/6 


COS  0)(u>  •  u>)  , 


2  -s-1/6 


=  Cj  +  2. A3  ^1  -  -j  cos2  0^  (w  •  w) 


substituting  (66)  into  (65), 


(w,  Aw,  Ay)  = 


- 5 — 

2*  1,2  (  2 ,  \  . 

j -  exp  <-k  y  [U2  ~  PL/uJ  Ay  • 

<k  Vo 

f  2  ^  2  1 

exp  < -k  Y^V  Aw  •  Aw  f  exp  <  -k  c  y  > 
i  2o  j  [  oo) 


(67) 


(68) 


(69) 


From  (69)  we  see  that  Mq(w,  Aw,  Ay)  behaves  as  a  Gaussian 
function  in  both  the  Ay  and  Aw  directions.  We  may  take  advantage 
of  this  by  defining 
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-1 


,  2 

k  Vo 


-1 


Equation  (69)  becomes 


Mo(u),  Aw,  Ay) 


2tt 


kVo 


,21 


exp  |~Ay  •  Ay/Ay^j>  exp  < -Au>  •  Aa>/Aco^ 


2 1 


•  exp 


/  ,2  ) 

s  -k  c  p 

l  0  °J 


(70) 


Equation  (70)  involves  only  a  quadratic  term  in  Ay  and  hence 
the  Fourier  transform  implied  by  (37)  can  be  calculated  in  closed 
form.  The  result  is 


M  (w,  Am,  Aw' )  =  //  M  (w,  Au),  Ay)  exp  [jkAw’  •  Ay]  d  Ay 


2* 


Ay 


Jc  j  2 

- - exp  ^  -Au  ♦  Au/Ati)  :■ 

u  i  CJ 


k  Vo 


I  2  2  I 

v  — k  Av  Adi  *  •  Am1  /L> 


l 


f  .  2 


exp  -k  Ay^Aoi*  •  Au’/Ay  exp  <^-k  • 


o  o 


(71) 


Defining 


<*>;  5  2/kAyc, 


(72) 
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I 


l 


we  have 


Mq(w,  Aw,  Aw')  =  2n‘l 


Ayc 

— - exp  -Aw 

k  Vo 


Aw/  Au> 


2 

c ) 


exp 


-Aw* 


Aw,/(Aw^)2| 


exp 


f  2 

■  -k  c  |i 
(.  o  o 


}• 


(73) 


Hence,  we  have  found  approximate  analytical  expressions  for 

Mo(w,  Aw,  Ay)  and  Mo(w,  Aw,  Aw1). 

We  now  consider  the  problem  of  selecting  the  coefficients  c^ 

and  e  .  These  coefficients  are  selected  so  that 
o 


d(p)  “  cJp 


P  +  c 

o 


(51) 


when  p  »  w^.  The  coefficients  cQ  and  c^  were  calculated  and  com¬ 
pared  using  a  variety  of  methods.  The  details  of  these  calcula¬ 
tions  are  relegated  to  Appendix  C.  The  following  values  for  Cj  and 
cQ  were  found  to  be  reasonable  choices: 


c 


1 


(74) 


c  =  0. 

o 


(75) 


Substituting  these  values  into  our  previous  results,  we  obtain 


M  (w,  Aw,  Ay)  “ 
o 


■  —  exp  -Ay  •  Ay/Ay2*  exp  <j-Aw  •  Aw/  Aw2  j' 


k  Vo 


(76) 


1 


30  - 


and 


with 


M  (w,  Aw,  Aw' )  »  2ti‘ 

o 


Ay 


c  ....  I  .  .  ,.  2  l 


r  v 


exp  -  -Aw  •  Aw/ Aw 


c  J 


(  9 1 

*  exp  ^-Aw'  •  Aw'/(Aw^)  j’  (77) 


i  ’ 12  (A„f  - 9-73  (» -  j  co“2  8) 


2  0\,“  -v-1/6 

cos  6](w  •  w)  , 


(78) 


Y2  =  3  (k2wo^  +  2. A3  (l  -  ~  cos2  0^  (w 


-v-1/6 

“> 


Aw  =  -*? 

c  1.2. 


k  Vo 


(79) 


(80) 


Ay, 


(81) 


and 


Aw' 

c 


■  2  |h(v  "i'O- 


(82) 


The  validity  of  the  approximations  used  to  derive  (76)  and 
(77)  were  tested  by  comparing  the  value  of  M£)(a)»  Aw,  Ay)  predicted 
by  (76)  with  the  value  computed  by  numerical  integration  of  (30). 
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The  two  results  were  compared  for  a  family  of  values  of 

w,  Aw,  and  Ay.  The  details  of  the  numerical  analysis  and  the  tech¬ 
nical  aspects  of  the  numerical  integration  are  included  in  section 
6  and  Appendix  D  of  this  report.  However,  initial  numerical  re¬ 
sults  showed  that  (76)  closely  approximated  the  value  obtained 
through  numerical  integration  of  (30)  for  large  |w|  (|w|  >  10w^). 
As  |w|  was  decreased  (76)  became  an  increasingly  poor  model  for 
Mo(w,  Aw,  Ay).  The  model  completely  broke  down  as  |w|  +  w^ .  The 
approximations  used  to  derive  (76)  assumed  that  |w|  >>  w^ ;  hence, 
it  is  not  surprising  that  the  expression  for  M^tw,  Aw,  Ay)  ceases 

to  be  accurate  for  |w|  *>  w  .  We  will  analyze  this  problem  in  more 

c 

detail  to  find  the  lower  bound  on  |wt  for  which  our  model  is  valid. 
This  lower  bound  will  quantify  the  previous  somewhat  nebulous  re¬ 
quirement  that  |w|  >>  w  . 

c 

We  return  to  (34).  Aw,  Ay)  was  found  by  integrating 

M'(w,  Aw,  A0,  Ay)  with  respect  to  A0  over  the  infinite  A0  plane. 
We  should  note  here  that  the  integration  over  the  aperture  plane  is 
actually  limited  by  the  size  of  the  aperture.  Thus,  physically  the 
integral  always  converges.  However,  in  deriving  (35)  we  argued 
that  the  effects  of  the  finite  apertt  ce  could  be  separated  from  the 
integration  over  A0  by  sampling  P'  (w,  Aw,  0,  A0)  at  A0  =  -Aw/2. 
We  now  explore  the  consequences  of  this  approximation  and,  in  par¬ 
ticular  we  determine  the  range  of  w  for  which  the  approximation  is 
reasonable.  For  the  integral  in  (34)  to  exist,  M'(w,  Aw,  A0,  Ay) 
must  vanish  for  sufficiently  large  |A0|.  We  shall  explore 
the  behavior  of  M'(w,  Aw,  A0,  Ay)  for  very  large  A0.  From  (30) 
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r 

:  i 


it  is  apparent  that  for  very  large  AB,  the  terms  d^(«), 

d^(*),  and  d^(*)  approximately  cancel  each  other  and  we  have 


M'(w,  Aoo,  AB,  Ay)  *  exp  <- 


f  « 


/  C  (h)  [d(u)  +  Aw)  +  d((o  -  Aw)]  dh' 

o  n  i 


|AB|  very  large. 


(83) 


We  are  not  interested  here  in  the  Aw  dependence  so  we  set  Aw  =  0. 
Substituting  (23)  into  (83), 

f  ,  H  . 

M'(w,  Aw,  AB,  Ay)  -  exp  --2.92  k  /  C  (h)  dh  |w|' 

i  o  n 

!  2  -  S/3  ^ 

=  exp  '-2.92  k  vj  \u\  1  >  .  (84) 

t  o  J 

From  (84)  we  observe  that  if  |w!  is  not  large  enough ,  then 
M'(w,  Aw,  AS,  Ay)  will  not  effectively  vanish  for  large  !  A3 1  and 
the  integration  over  the  AS  plane  is  not  defined.  We  hasten  to  add 
that  (84)  indicates  that  M'(w,  Aw,  AB,  Ay)  does  not  strictly  vanish 
for  large  |ABl  but  rather  approaches  some  constant  value.  We  must 
quantify  what  we  mean  by  effectively  vanishing.  Numerical  analysis 
has  shown  good  correlation  between  (76)  and  the  numerical  integra¬ 
tion  of  (30)  if  |w|  is  large  enough  to  insure  that  the  asymptotic 
value  predicted  in  (84)  is  on  the  order  of  10  ^  or  smaller.  We  can 
now  compute  the  lower  bound  on  |w|  subject  to  this  condition.  The 
lower  bound  will  be  denoted  by  w^,  and  is  calculated  from  the 
condition 
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exp 


-2.92  k  u 


5/3 


cl 1 


=  10 


-6 


Solving  for  to  , , 
cl 


we  find  that 


o> 


13.8 


\3/5 


cl 


,2.92  k 


3  to 


(85) 


We  have  established  to  ,  as  a  lower  bound  on  I  to [ .  Our 

cl 

model  for  M  (to,  Am,  Ay)  is  valid  for  |to|  >  to  ..  We  can  also  inter- 
o  cl 

pret  physically  the  behavior  of  M’(io,  Ato,  Af5,  Ay)  in  the 

region  to  <  { to (  <  to  In  section  4  we  noted  that  for  spatial 

c  cl 

-  2  2  2 

frequencies  satisfying  |  to  |  +  1  Ato)  <  u>c,  the  entire  aperture 
appears  to  be  spatially  coherent.  Furthermore,  for  spatial  frequen¬ 
cies  satisfying  I  Am|  <  uj^  and  l to l  »  co^ ,  the  aperture  is  spatial¬ 
ly  coherent  in  regions  the  size  of  to  (for  source  points  separated 

by  Ay  <  Ay  ).  The  region  to  <  |to|  <  to  ,  is  in  the  transition 
c  c  cl 

zone  between  the  regions  of  complete  and  partial  coherence. 

There  is  one  remaining  issue  that  must  be  resolved.  In  the 
course  of  deriving  (76)  and  (77),  we  assumed  that  l  to l  >  | h Ay |  , 
knowing  that  this  assumption  was  suspect.  The  preliminary  numeri¬ 
cal  analysis  mentioned  above  has  shown  additional  discrepancies 
between  (76)  and  the  values  calculated  from  numerical  integration 
of  (30)  for  values  of  |to|  larger  than  w  The  region  in  which  the 


discrepancy  was 

significant 

was 

generally  confined 

to 

to  .  <  j  to  |  <  5to  .  . 
cl  cl 

For 

larger  v  . ' 

ues  < 

of  |<o|,  (76)  was  valid. 

We 

should  point  out 

that 

the  width 

of 

this  troublesome  region 

is 
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highly  dependent  upon  the  model  used  for  the  refractive  index 

2  2 
structure  constant  C  (h).  The  width  is  largest  for  C  (h)  profiles 

n  &  n 

with  very  weak  turbulence  at  higher  altitudes. 

The  discrepancy  is  in  the  Ay  plane  only.  For  values  of  ui  in 

the  troublesome  region,  M  (w,  Ato,  Ay)  no  longer  appears  to  be  a 

o 

Gaussian  in  Ay  but  rather  the  sura  of  a  Gaussian  and  a  constant. 
The  value  of  the  constant  decreases  with  increasing  |a>!.  The 
problem  is  due  to  the  fact  that  under  these  conditions  |w|  does  not 
dominate  |hAy|  and  the  truncated  binomial  expansions  used  in  the 
derivation  of  (76)  are  not  accurate. 

To  quantify  the  effects  we  have  observed,  we  seek  an  alter¬ 
native  approximation  to  M'(u),  Au>,  A8,  Ay)  that  is  valid  for  large 
1 hAy |  .  Unfortunately,  we  have  not  been  successful  in  determining 
an  approximation  that  preserves  integrability  with  respect  to  h. 

Hence,  we  must  sacrifice  generality  here  and  assume  some  model  for 
2 

the  C  (h)  profile.  To  expedite  the  integration  over  h,  we  will  use 
n 


0  h1  h 

2 

Fig.  2.  The  double  impulse  model  for  C^Ch). 
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We  model  C  (h)  by 
n 


C2(h) 

n 


-  aQ6(h)  +  a^6|_h  -  h^), 


a  >  a,  . 
o  1 


The  impulse  at  h  =  0  represents  turbulence  near  the  ground  while 
the  impulse  at  h  =  hj  represents  turbulence  in  the  upper  atmo¬ 
sphere.  The  height  hj  of  the  second  impulse  is  in  general  vari¬ 
able,  but  in  this  report  we  will  associate  with  the  "tropopausal 
bump"  sometimes  observed  at  a  height  of  10^  m  [8],  VJe  require 

that  a  >  a,  to  satisfy  (12). 

o  1 

The  double  impulse  model  allows  the  trivial  evaluation  of  the 
integration  with  respect  to  h.  Substituting  the  double  impulse 
model  into  (10)  and  integrating  over  h  yields 


M'(m,  Am,  AS,  Ay)  = 


/  8^(0+  Am)  +  d  (m  -  Am) 

k2 

exp  -~  al  +  d3(2A3  +  2Am)  +  d^(2AB) 


-  d^(2AB  +  Am  -  u)  -  d^(2AB  +  Am  +  m). 


d7(m  +  Am)  +  dQ(m  -  Am) 

/  O 


+  I  +  dg(2A8  +  2h^Ay  +  2Am)  +  d^p(2AB  +  2h^  Ay)  J 
\  ^1  j (2AB+2h^ Ay+Am-U')  -  d j2(2A8+2h^Ay+Am+m)/ 


J  J 

(87) 
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Once  again  we  will  approximate  each  of  the  d(»)  functions  by  qua¬ 
dratics.  However,  here  we  are  interested  in  large  values  of 
Ih^Ayl,  so  we  will  assume  that  Ih^Ay]  dominates  the  other  terms  in 
the  arguments  of  d^(»)  to  d  (•).  Previously,  we  have  shown  that 

d(p)  *  c ,  p  •  p  +  c  fp  •>  u  )  (51) 

1  o  c 

and 

d(u  +  6)  »  2.92  (w  •  oi)5/6  1  +  ~~"r  +  ~~*'V ( f- 

3ui  •  w  a)  •  u>  ' 

Similarly, 

/  9  \5^6  5h  5  •  Ay  A  .  6 

d(2h.  Ay  +  6)  *  2.92  (  Ah- Ay  •  Ay  ]  1  + - 2 -  + - 5 - 

'  12h^Ay  •  Ay  Ah^Ay  •  Ay 

*  (|  ‘  if  C°sZ  *)  *  (88) 

where  \|j  is  the  angle  between  Ay  and  5.  We  are  interested  here  in 
the  qualitative  behavior  of  M1  (c*>,  Au),  AB,  Ay)  as  1  h Ay |  becomes 
very  large.  To  simplify  the  analysis,  we  will  neglect  the  angular 
dependence  and  set  cosz  =  1.  We  approximate  d^(*)  and  d^(»)  by 
(51),  d^  ( • ) ,  d^(* ) ,  d^(* ) ,  d&(*),  and  dg(*^  by  (46), 

and  d^( • ) ,  d1Q(*),  d^f*),  and  d12(*)  by  (88).  In  Appendix  C, 
cQ  is  shown  to  be  nearly  zero.  Making  these  substitutions  in  (87), 


we  obtain 


M’(a),  Am,  A6,  Ay) 


/  4c.  -  6.49(wid)  1/6,a  A$‘A£5+  4c  -6 .49(id*w) 
/  l  1  I  0  1 


I  a  Aw* 
I  o 


As\ 


!  ,2 

exp  -k 


+  2c,  a  +  1.62(a)  •  a))  | 

1  o  lj 


Aio  •  Au) 


+  |2 . 92 (d)»d))  -  1 . 62  /4h2 Ay  •  Ay  \  ( d)*d)-Ad)*  Aid) 


(89) 

This  equation  is  quadratic  in  Af3.  Completing  the  square  and  inte¬ 
grating  over  the  A8  plane  yields 


M^(d) ,  Ad),  Ay) 


~ —  exp  -k2  3fk2p  )  a  +  1.62(d)  •  u) 

!  LlVo/o  o. 


k2y'o 


-n-1/6  [  '  .  \ 

d))  p_  |  Ad)  •  Ad)  r 


[of  __  __  C /A  {  7  1/6  __  __  -i'l 

j^-k  1^2 .92(d)*w)  -  1.62  ^4h^Ay*  AyJ  (d)*d>-AH)»  Ad))  > 


exp 


J 

(90) 


where  Y ’  is  given  by: 


1/5 


Y'  =  12  (k2p  )  -  6.49  (id  «  i)"1/6. 


(91) 


The  Ay  dependence  is  no  longer  Gaussian,  and  we  emphasize  that  this 
approximation  is  valid  only  for  very  large  |  h  Ay  |  .  In  particular, 
setting  Ad)  =  0  and  taking  the  limit  as  |h^Ay|  ♦  <*>, 


lim  MQ(d),  0,  Ay) 

|h1Ay I +« 


2ir 


,  2  , 
k  Y  a 


exp 


f  2 

j-2.92  k  ax 


|Sl5/3} 


(92) 
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Equation  (76)  predicts  that  H^(u,  0,  Ay)  +0  as  | h ^ Ay  1 
However,  (91)  predicts  that  0>  Ay)  g(w,  ,  h^  )  *  0 

as  |h^  Ay|  •»  ®.  We  also  note  that  as  |id|  increases,  the  value  of  g 
decreases.  To  bound  the  region  where  this  effect  is  significant, 
we  propose  to  neglect  g  when  |i»)|  >  us^-  We  define  u  ^  as 


exp  j-2.92  k  l“c2i 


5/3 


exp 


{-2} 


(93) 


Solving  for  yields 


\3/5 


u  = 
c2 


.2.92  k  a. 


(94) 


It  is  more  convenient  to  express  to  as  a  multiple  of  w  .  Using 

C2  c 

(27),  we  obtain 


(93) 


Equation  (95)  bounds  the  region  where  the  noted  phenomenon  is  sig¬ 
nificant.  The  effect  is  significant  for  spatial  frequencies  id 

such  that  u  <  |u>|  <  u>  .  We  note  that 
cl  c2 


u  ,  “  3  w 
cl  c 


(96) 


u>  _  *  Ku> 
c2  c 


(97) 


where  K  depends  upon  the  ratio  of  surface-to-high-altitude  turbo- 
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lence.  Some  values  of  K  are  tabulated  in  Table  1. 


Table  1.  Representative  values  of  K  as  a  function  of 

the  rat io  a  /a, . 

o  1 


a  /a 
o'  1 

K 

100 

15.9 

10 

4.2 

8 

3.7 

4 

2.3 

2 

1.5 

- - * 

We  now  make  a  key  observation.  If  a  /a,  <  5.2  (i.e.,  if  a, 

o  1  1 

is  at  least  20  percent  of  a  ),  then  ui  „  <  o>  ,  and  (76)  is  valid  for 

o  c2  cl 

all  | cu |  >  oi  ,  .  On  the  other  hand,  if  a  /a  >  5.2,  then  there  will 
cl  o  1 

be  a  nonvanishing  region  in  which  the  Gaussian  model  in  Ay  must  be 

replaced  by  a  Gaussian  plus  constant  model.  However,  the  width  of 

this  region  does  not  become  significant  until  o  /a,  >  10. 

o  1 

At  this  point  it  is  useful  to  partition  the  w  plane  into  four 
regions : 

a.  1  oi |  <  oi 

c 


cu  <  |  oi|  <  oi  , 
c  cl 


c.  0)cl  <  | Oil  <  0)c2 


-  A0  - 


d. 


“c2  <  l“l 


Under  the  assumptions  |Aw|  <  |  2A6  +  Aw|  <  w^  and  [  Ay  [  <  Ay^ , 

in  region  (a)  the  entire  aperture  appears  coherent.  Region  (b)  is 

the  transition  zone  from  complete  to  partial  coherence.  In  region 

(c),  ( to ,  Aw,  Ay)  is  best  modeled  by  the  product  of  a  Gaussian  in 

Aw  and  a  Gaussian  plus  constant  in  Ay.  In  region  (d) 

M^iw,  Aw,  Ay)  may  be  modeled  by  the  product  of  a  Gaussian  in  Aw 

and  a  Gaussian  in  Ay.  We  have  shown  that  region  (c)  vanishes  for 

moderate  to  strong  upper  atmospheric  turbulence. 

Equations  (89)  and  (90)  are  strictly  valid  only  for  the 
2 

double-impulse  model  of  C^(h).  However,  the  analysis  should  extend 
2 

to  more  general  C^(h)  profiles.  The  key  result  is  contained  in 

(92).  We  conclude  for  the  general  case  that  if  some  appropriate 

measure  of  the  turbulence  in  the  upper  atmosphere  is  much  less  than 

the  corresponding  measure  of  turbulence  near  the  surface  (i.e., 

less  than  10  percent),  then  there  will  exist  a  nonvanishing  region 

(c)  w  <  |w|  <  w  where  our  models  for  M  (w.  Aw,  Ay)  and 
c  1  c  i  o 

Mq(w,  Aw,  Aw')  given  by  (76)  and  (77)  will  not  be  valid.  For  the 

2 

double-impulse  model  of  C^Oi),  we  can  predict  the  size  of  the 
region  and  the  form  of  the  deviation  from  (76)  and  (77).  In  the 
more  general  case,  the  exact  form  of  the  deviation  cannot  be  pre¬ 
dicted  by  a  closed  form  expression.  We  note  that,  physically,  the 
phenomenon  we  have  described  is  the  transition  of  the  imaging 
problem  from  the  isoplanatic  case  to  the  nonisoplanat ic  case. 

One  reasonable  technique  for  measuring  the  relative  amounts 
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of  upper  atmospheric  and  surface  turbulence  for  an  arbitrary  sample 
2 

function  of  C  (h)  is  to  compute  the  moments  p  ,  P.  ,  and  p  of  the 
n  o  i  i 

sample  function  and  then  convert  these  values  into  an  equivalent 
double-impulse  model.  Equation  (95)  can  then  be  applied  to  deter¬ 
mine  the  approximate  value  of  F°r  example,  consider  Huf- 

2 

nagel's  simple  model  [9]  for  C^ih): 


C2(h) 

n 


h  <  20,000 
h  >  20,000 


(98) 


The  moments  of  this  profile  can  be  easily  computed: 

p  -  1.64  ♦  10"12, 
o 

=  3.00  •  10-9, 

P2  =  3.00  •  10-5. 

An  equivalent  double-impulse  profile  would  have  parameters 

o  =  1.34  •  10"12, 
o 

a1  -  3.00  •  10”13, 

4 

ht  =  10  . 
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Computing  u>  ,  we  find  that 
cl 


c  2 


2.8  w  . 

c 


We  conclude  that  region  (c)  vanishes  for  this  profile  and  (76)  and 
(77)  are  valid  approximations  for  all  |m|  >  3^. 

In  general,  region  (c)  is  significant  only  for  atmospheric 

models  with  extremely  weak  high  altitude  turbulence.  From  the 

above  analysis  of  Hufnagel's  simple  model  and  other  published  pro- 
2 

files  of  (^(h),  it  appears  reasonable  to  conclude  that  models  with 
such  extremely  weak  high  altitude  turbulence  are  not  often  encoun¬ 
tered  in  practice. 

Returning  to  the  questions  posed  before  embarking  upon  this 
series  of  derivations,  we  find  that  simple  analytical  expressions 
can  be  found  that  approximate  Mo(w,  Au>,  Ay)  and  Am,  Aw'). 

Our  models  for  these  two  functions  are  given  by  (76)  and  (77).  We 
find  that  C^(h)  influences  M^Cu),  Am,  Ay)  and  Am,  Am')  through 

the  zero,  first,  and  second  moments  (pq,  ,  and  These  are 

the  salient  parameters  of  the  refractive  index  structure  constant, 
and  aside  from  condition  (12)  the  general  shape  of  the  sample  func¬ 
tion  is  unimportant.  We  will  discuss  the  impact  of  the  smoothing 
function  Mo(w,  Aw,  Au>')  on  the  imaging  statistics  in  section  8  of 
this  report. 


6.  Numerical  Verification  of  the  M„  Model 

In  this  section  of  the  report,  we  will  summarize  the  results 
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of  numerical  analysis  performed  to  test  the  validity  of  the  closed 
form  expression  for  MQ.  The  integration  over  A8  called  for  in  (34) 
was  computed  numerically  and  compared  with  the  value  predicted  by 
our  model.  Only  one-dimensional  problems  were  considered  to  keep 
the  numerical  routines  tractable.  The  details  of  the  numerical 
integration  algorithm  are  discussed  in  Appendix  D.  The  one-dimen¬ 
sional  equations  for  M0(w,  Am,  Ay)  are  given  in  Appendix  E. 

The  calculation  of  Mq(w,  Aw,  Ay)  requires  two  integrations: 
the  first  with  respect  to  h,  and  the  second  with  respect  to  AB. 

Initially,  the  numerical  complexity  of  the  integration  routine  was 

2 

simplified  by  presuming  the  double  impulse  model  for  C  (h)  (Fig. 

n 

2).  This  supposition  makes  the  integration  over  h  trivial. 

Three  double-impulse  models  were  considered.  They  differ  in 

the  ratio  of  upper-to-lower  atmospheric  turbulence  (i.  e., 
2  2 

C  fh. )/C  (0)).  The  model  parameters  are  summarized  in  Table  2. 
nv  1 '  n 

Table  2.  Parameters  of  the  three  double-impulse  models. 


A 

C  (h)  model  \ 
n 

a 

o 

— 

al ' 

hl 

i% 

10-12 

10'14 

10* 

10% 

10"12 

10-13 

104 

25% 

10-12 

2.5  •  10~13 

1 

O 

We  refer  to  the  three  profiles  as  the  1  percent,  10  percent,  and  25 
percent  models,  where  the  percentage  is  the  ratio  of  upper-to-lower 
atmospheric  turbulence  ((o^/a^)  •  100).  In  each  case  we  set  the 
height  of  the  second  impulse,  hj ,  at  10,000  meters.  The  various 


atmospheric  cutoff  frequencies  f to  ,  w  ,  ,  w  1  for  the  three  C  (h) 

n  1  c  cl  c2  ‘  n 

profiles  are  tabulated  in  Table  3. 

Table  3.  Cutoff  frequencies  for  the  three  double-impulse  models. 


C2(h)  Model 
n 

D 

C 

W  1 
cl 

V 

1% 

0.050 

0.150 

0.797 

10% 

0.047 

0.141 

0.198 

25% 

0.044 

0.132 

0.116 

All  three  models  have  comparable  values  of  w  and  w  .  However, 

c  c  1 

the  value  of  w^  varies  considerably  among  the  three  models.  This 
is  due  to  the  wide  variation  in  the  relative  amounts  of  upper- 
atmospheric  turbulence  in  the  three  models.  The  1  percent  model 
implies  extremely  weak  upper-atmospheric  turbulence  and  hence  the 
value  of  <0^2  is  large.  The  10  and  25  percent  models  imply  moderate 
and  moderate-to-strong  upper-atmc spheric  turbulence  and  have  corre¬ 
spondingly  smaller  values  of  w^.  Indeed,  for  the  25  percent 
model,  region  c  has  vanished  (section  5). 

The  values  of  Mq(w,  Aw,  Ay)  computed  from  (E2)  and  by  the 
numerical  integration  of  (30)  were  compared  for  many  (350)  combina¬ 
tions  of  w,  Aw,  Ay.  The  correlation  was  found  to  be  very  good, 
with  errors  generally  less  than  10  percent  until  M^iw,  Aw,  Ay) 
approached  zero. 

Some  typical  results  are  plotted  in  Figs.  3-5.  In  each  plot 
cross  sections  in  the  Aw  and  the  Ay  directions  are  shown,  with  the 
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other  parameters  fixed.  The  numerical  integration  data  are  plotted 
in  a  dashed  line,  while  the  values  predicted  by  ( E2 )  are  plotted  in 
a  solid  line. 

Figures  3a,b-5a,b  illustrate  the  roll-off  of  M^(w,  Aw,  Ay) 
in  the  Aw  direction  with  Ay  =  10  ^  rads.  Since  Ay  <<  Ay^ ,  there  is 
effectively  no  attenuation  of  Mo(w,  Aw,  Ay)  due  to  the  separation 
Ay  of  the  image  points.  Similarly,  Figs.  3c,d-5c,d  illustrate  the 

_  _3 

roll-off  of  M^(w,  Aw,  Ay)  in  the  Ay  direction  with  Aw  =  10 

meters.  Since  Aw  <<  w  ,  there  is  effectively  no  attenuation  of 

c 

M^(w,  Aw,  Ay)  due  to  the  difference  in  aperture  spatial  frequen¬ 
cies  Aw.  Figures  6a-6b  are  analogous  to  5a  and  5b  except  that  a 
much  larger  value  of  Ay  was  selected  (Ay  =  3  •  10  to  illustrate 
the  cross  section  in  Aw  when  Ay  is  sufficiently  large  to  contribute 
to  the  attenuation.  Also,  Figs.  6c  and  6d  are  analogous  to  5c  and 
5d  except  for  the  selection  of  a  larger  value  of 

Aw  (Aw  *  3  •  10  ),  to  illustrate  the  cross  section  in  Ay  when  Aw 

is  sufficiently  large  to  contribute  to  the  attenuation. 

Figures  7a  and  7b  demonstrate  the  rapid  degradation  in  the 

model  for  M  (w  Aw,  Ay)  (E2)  for  w  <  w  . .  From  Table  3,  for  the  25 
o  cl  ’ 

percent  profile,  w^  =  0.116  meters.  Figures  7a  and  7b  are  plots 
of  the  Aw  and  Ay  cross  section,  with  w  =  0.100  meters.  For  this 
value  of  w  slightly  less  than  w^,  the  degradation  in  our  model  is 
already  apparent.  Other  numerical  results  indicate  that  the  break¬ 
down  of  the  model  is  rapid  and  severe  for  w  <  w^. 

There  is  strong  motivation  to  further  verify  the  model  for 
Mq(w,  Aw,  Ay)  by  generalizing  the  numerical  comparisons  made  above 
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2 

to  more  general  sample  functions  of  C^(h).  Unfortunately,  the 
computational  complexity  of  the  "nested"  integration  implied  in 
(30)  can  rapidly  become  overwhe lming  even  for  reasonably  "well- 
behaved"  profiles.  Rather  than  tackle  this  more  comprehensive 
problem,  we  have  tested  a  sampled  version  of  Hufnagel's  simple 
model  : 


C2(h) 

n 


h  <  20,000 


|  0  h  >  20,000 

This  profile  was  sampled  at  1  km  intervals  resulting 
pie  impulse  model  shown  in  Fig.  8  with 


in 


(99) 


the  multi- 


a 

o 


500 

»  i 

0 


1.5 


10 


-13 


dh , 


(2i+l)500 


ai  =  J 


1.5 


10 


■13 


dh 


(2i-l)500 


1  <  i  <  19, 


20 


20,000 

=  i 

19,500 


1.5 


10 


■13 


dh . 


The  corresponding  Aw  and  Ay  cross-sectional  plots  generated  by  the 
numerical  integration  of  (34)  with  this  profile  and  the  values 
predicted  by  (El)  are  shown  in  Figs.  9a-9d.  This  model  has  a  sig¬ 
nificant  level  of  upper-atmospheric  turbulence,  and  so  (El)  was 
used  rather  than  (E2).  Once  again,  the  dashed  line  represents  data 
calculated  by  numerical  integration  while  the  solid  line  represents 
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Fig.  3.  A  comparison  of  (E2)  with  the  numerical  integration  of  (34)  for 
the  1  percent  Cf^(h)  profile  (solid  line  -  (E2);  dashed  line  - 
numerical  integration). 

(a)  u  =  0.2,  Ay  =  10“5 

(b)  J.  =  1,  Ay  =  10"' 

(c)  "  =  0.2,  A_j  =  10'^ 

(d)  "  =  1  ,  Ai,:  =  10"-5 


» 


Fig.  4.  A  comparison  of  (E2)  with  the  numerical  integration  of 
the  10  percent  Cp(h)  profile  (solid  line  -  (E2);  dashed 
numerical  integration). 

(a)  H)  =  0.2;  Ay  =  10  7 

(b)  u  =  1,  Ay  =  10"o 

(c)  w  =  0.2;  Ago  =  10  ? 

(d)  w  =  l;  Au  =  10'J 


(34)  for 
line- 


V^v  .  .  . 


L)i  .  r  a  r 
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Fig.  5.  A  comparison  of  (E2)  with  the  numerical  integration  of  (34)  for 
the  25  percent  Cfi(h)  profile  (solid  line  -  (E2);  dashed  line  - 
numerical  integration. 

,-7 


(a) 

y  = 

0.2, 

Ay  = 

10 

(b) 

U)  - 

1, 

Ay  = 

10 

(c) 

y  ~ 

0.2, 

Aoo  = 

10 

(d) 

U)  = 

1, 

Aui  = 

10 
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Fig.  6.  A  comparison  of  (E2)  with  the  numerical  integration  of  (34)  for 
the  25  percent  c£(h)  profile  (solid  line  -  (E2);  dashed  line  - 
numerical  integration. 


(a)  u  =  0.2,  Ay  =  3 

(b)  uj  -  1 ,  Ay  =  3 

(c)  u  =  0.2,  Aw  =  3 

( d )  w  =  1 ,  Aw  =  3 
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Fig.  7.  A  comparison  of  (E2)  with  the  numerical  integration  of  (34)  for 
the  25  percent  C£(h)  profile  with  I  <  wc-|  (solid  line  -  (E2); 
dashed  line  -  numerical  integration). 

(a )  u  =  0.1 ,  Ay  =  1 0  \ 

(b)  w  =  0.1 ,  Aui  =  10" 
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ig.  9.  A  comparison  of  (El)  with  the  numerical  integration  of  (34)  for 
the  sampled  Hufnagel  model  C^(h)  profile  (solid  line  -  (El); 
dashed  line  -  numerical  integration. 

(a)  w  =  0.2,  Ay  =  10_^ 


our  model. 


In  conclusion,  we  have  shown  the  closed  form  expression  (E2) 
for  Mq(w,  Aw,  Ay)  to  be  a  reasonable  approximation  by  comparing 
(E2)  to  data  obtained  by  numerical  integration  of  (34).  In  partic¬ 
ular,  our  numerical  experiments  verified  the  following  phenomena: 

1.  The  behavior  of  Mq(w,  Am,  Ay)  in  the  Aw  direction  is 

modeled  well  for  all  w  >  w  ,  by  a  Gaussian  curve. 

cl 

2.  The  behavior  of  Mq(w,  Aw,  Ay)  in  the  Ay  direction  is 

modeled  well  for  w,  <  w  <  w  „  by  a  Gaussian  curve  plus  a 
constant.  For  w  >  w^)  t^ie  behavior  of  Mo(w,  Aw,  Ay)  may 
be  modeled  by  a  Gaussian  curve  in  Ay.  The  approximation 

improves  for  increasing  w. 

2  ■  .  , 

3.  For  C^(h)  profiles  with  significant  upper-atmospheric 

turbulence,  w  „  <  w  ,  and  we  may  neglect  the  phenomenon 
cZ  c  1  ' 

described  in  (2). 

4.  For  w  <  w  ,,  the  model  ceases  to  be  valid.  The  breakdown 

cl 

of  the  model  in  this  region  is  rapid  and  severe. 


7.  Comparison  with  Fried 

In  this  section  of  the  report,  we  will  compare  some  of  the 
implications  of  our  model  for  Mq(w,  Aw,  Ay)  with  similar  conclu¬ 
sions  made  by  Fried,  et  al.  [10,  11].  We  can  associate  the  param¬ 
eter  Ayc  of  the  model  with  the  size  of  the  isoplanatic  patch. 
Fried  has  proposed  the  following  expression  for  Ay^ : 


Ay«  = 


2.91  k2  / 


h^3  C2(h)  dh 
n  l 


-3/5 


(100) 
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We  have  compared  the  isoplanatic  patch  sizes  predicted  by  Fried 

2 

(100)  and  by  oar  model  (81)  for  the  C  (h)  profiles  discussed  in  the 

n 

previous  section. 


Ay,  = 


Jk2 

l 


P?/U  ) 

1  o  'J 


-1/2 


(81) 


with 


Yj  =  12  (k2Po)  -  9.73  (l  "  3  cos2  9^  (u  •  u)  1/6.  (78) 

Since  our  model  includes  some  weak  dependence  on  i  to  1  and  the  angle 
between  u  and  Aio,  we  will  set  |oi|  =  1  and  0  =  0  for  the  purpose  of 
comparison.  The  results  are  tabulated  in  Table  4. 


Table  4.  A  comparison  of  the  isoplanatic  patch  size 
predicted  by  (81)  and  by  Fried. 


- - — - - - 

C^(h)  Profile 
n 

Eq.  81 

Fried 

1%  model 

Ay  =  2.06  •  10-5 

1  c 

Ay*  =  5.27  •  10~5 
c 

10%  model 

6.74  ♦  10~6 

1.32  •  10~5 

25%  model 

4.47  •  10~6 

7.64  •  10"6 

Sampled  Hufnagel  model 

_  . . 

3.98  •  10“6 

7.03  •  10~6 

We  find  very  reasonable  correlation  with  differences  on  the  order 
of  a  factor  of  two. 
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8.  Effects  on  Imaging  Equations 

We  will  now  consider  the  effects  of  M  (w,  Aw,  Aw')  on  the 

o 

imaging  statistics.  We  assume  M^(w,  Aw,  Aw')  to  be  of  the  form 
given  in  (77): 


M  (w,  Aw,  Aw')  »  V  exp  {-Aw* 


Aw'  /  (Aw^_  } 


(101) 


where 


V 


2* 


,  2 

k  Vo 


2 

exp  {-Aw  •  Aw/(Aw^)  } 


(102) 


Substituting  this  expression  for  ^(w,  Aw,  Aw')  into  (39),  we 
deduce  the  following  expression  for  the  Fourier  transform  of  the 
second-order  statistics: 


wj  =  4|At 


4  (2it ) 


(k) 


V  F. (w, 
4 


•  -  / 

Aw)  JJ  0  (  w1  - 


:xp  {-Aw’  •  Aw'/(Aw')  }  d  Aw'. 


(103) 


In  the  following  analysis  we  shall  neglect  the  factor  preceding  the 
integral  in  (103).  This  factor  includes  scaling  constants  and 
gain-attenuation  factors  in  both  the  w  and  Aw  planes.  Our  primary 
interest  here  is  in  the  distortion  of  the  second-order  statistics 
in  the  Aw'  plane.  We  define 


» 
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(104) 


W'(wj.  “2)  s  II  0  ^  0  +*—- 


exp  { —  Aoj '  •  Aw'/{Aw')2}  d2Aw' 


To  provide  more  insight  into  the  smoothing  effects  of 
Aw,  Aw'),  we  shall  assume  the  image  to  be  Gaussian.  The 
Fourier  transform  of  a  Gaussian  is  still  a  Gaussian  so  we  may 
assume  that 


2 

0(w)  =  exp  {-w  *  w/.0  }  . 


(105) 


Substituting  (105)  into  (104),  we  obtain 


W  (<»!_,  u>2)  =  //  exp 


[  (w  -  Au'/2)  [  (u>2  +  Aw'/2) 

^ - - >  exp  < - j - 

0  J  l  0 


f\  o 

exp  {-Aa)'  •  Aa) ' / ( Aw^ )  }  d  Aw' 


(106) 


The  integral  may  be  computed  by  expanding  the  exponent  and  complet¬ 
ing  the  square  in  Aw' .  Neglecting  an  unimportant  scale  factor,  the 
integral  is 


|  2(Aw<|)2w  •  w  ] 

j(v  =  exp  ~2rT'"L;2 — exp  V 

1  2  [fi2[(Aw;)2  +  2021|  ( 


w  •  w 
1  1 


•  exp  •:  - 


w  •  w„ 

2  2  (. 


r  o2  )■ 


(107) 


We  observe  that  the  second  and  third  exponential  factors  in  (107) 
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are 


just  o(iOj  )  and  C^u^).  We  can  express  the  integral  j(w^,  w^ )  in 
terms  of  w  and  Aw  by  observing  that 


w1  •  w^  +  u>2  •  w2  = 


(<*>!  “  u2)  *  "  w2)  +  +  w2)  *  + 


2w  •  w  +  2 Aw  •  Aw. 


(108) 


Substituting  (108)  into  (107),  we  obtain 


J(w,  Aw)  =  exp 


2(Aw^)2  w 


'  2w  •  w  i 
?  exp  < - s —  > 


•  exp  s - 


1  r2[(aw*)2  +  2fi2]  j  l  n2  j 

\  2 Aw  •  AW' 

i"  *2  i  ' 


(109) 


Collecting  terms  in  w, 


J(w,  Aw)  =  exp 


2w  •  ui 


2n 


"\ 


) 


(Aw;)2+2n2]j  X?  [  r2  j 


2Aw  •  Aw 


(110) 


From  (110)  it  is  apparent  that  the  imaging  process  has  not  altered 
the  statistics  in  the  Aw  plane.  However,  the  imaging  process 
broadens  the  statistcs  in  the  w  plane.  Specifically,  the  width  of 
the  function  in  the  w  plane  is  increased  by  the  factor 

1  +  (Aw')2/2R2. 
c 
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Summary  and  Conclusions 


This  report  deals  with  the  problem  of  astronomical  imaging 
through  a  turbulent  atmosphere  under  conditions  that  force  us  to 
abandon  the  usual  assumption  of  isoplanic ity .  Integral  equations 
are  developed  which  relate  the  object  to  the  first  and  second  mo¬ 
ments  of  a  series  of  short-exposure  images.  These  equations  in¬ 
volve  the  structure  function  of  the  atmosphere,  which  is  shown  here 
to  be  reasonably  well  approximated  by  a  5/3-law  model  for  all 
values  of  its  argument. 

A  detailed  analysis  of  the  effects  of  nonisoplanatic  imaging 
potentially  involves  a  considerable  amount  of  numerical  computa¬ 
tion.  Much  of  this  computation  can  be  avoided,  and  at  the  same 
time  our  insight  into  the  process  can  be  enhanced  by  making  certain 
approximations  to  some  key  functions  which  describe  the  effects  of 
the  atmosphere.  Several  such  approximat ions  are  studied  in  this 
report,  and  their  validity  is  checked  against  results  obtained  by 
numerically  integrating  the  original  functions.  We  have  obtained 
some  useful  models  which  are  accurate  under  a  wide  range  of  reason¬ 
able  conditions  and  whose  use  greatly  simplifies  the  analysis  of 
nonisoplanatic  imaging. 
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APPENDIX  A 


In  this  appendix  we  consider  the  derivation  of  the  first-  and 
second-order  statistics  of  a  set  of  speckle  images  in  greater 
detail.  The  uncertainty  in  the  images  is  due  to  variations  in  the 
point-spread  function  s(x,  y).  To  obtain  expressions  for  the 
first-  and  second-order  statistics,  it  is  necessary  to  compute 
E{s(x,  y)}  and  e{s(x^,  y^ )  y2 ) }  -  from  (5)  we  find: 

oo  oo 

E{s(x,  y)}  =  |  A 1 2  jS  II  a(B  +  w)a(3)  exp[j0(B  +  to)  -  jB(B)] 

—  OO  — oo 

2  2 

•  exp[-jkto  •  (x  -  y)]  E{expt>K6  +  to,y)  +  t|>*(B,y)]}d  Bd  to 


where 


(Al) 


i(*,  •)  =  x(*»  ’)  +  j4>(*>  *).  (A2) 

We  assume  that  x(*j  *)  and  *)  ate  uncorrelated,  jointly  Gaus¬ 

sian  random  processes.  Furthermore,  wt-  assume  that 

e{x(*,  y)>  =  X  (A3) 


and 


E C 4> C  • ,  y )}  =  0. 


(AA) 
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For  notational  simplicity  we  define 


g(6,  gj;  y)  =  1^(1?  +  gj,  y)  +  ^*(8,  y). 


We  note  that  g(8,  gj;  y)  is  a  Gaussian  random  process. 

E  (exp[g]}  we  use  the  fact  that  for  a  Gaussian  random 

J  2 

with  mean  p  and  variance  0 


E{exp[x]}  =  exp 


Applying  (A6),  we  find  that 

E{exp[v|)(8  +  g>,  y)  +  ^*(8,  y)]}  =  exp  {E[g(8,  gj; 

+  E[(g(S,  gj;  y)  -  E(g(8,  gj;  y)))2]}. 


From  (A2-A5)  we  infer  that 


E{g(B,  g>;  y)}  =  2\(j 


and 


E{  [g(8,  gj;  y)  -  E(g(8,  gj;  y))]2}  =  e{[x(B  +  gj;  y)  + 

_  2 

+  j<}>(6  +  gj;  y)  -  j$(6,  y)  “  2x1  }  • 


( A5 ) 

To  compute 
variable  x 

(A6 ) 

y)l 

(  A7  ) 


(A8) 

X(8,  y) 

(A9) 


Since  x(*>  *)  and  $(*, 


)  are  uncorrelated  and  <}>  =  0,  we  can  elimi- 


nate  the  cross-product  terms  involving  both  x(*>  *)  and  (};(*,  •). 


Hence 

E{lg(B,  id;  y)  -  E(g(6,  id;  y ) )  1 2 }  =  e{(x(B  +  id,  y)  +  x(B,  y) 

-  2x  1  2 }  ~  e|[<>(3  +  id,  y)  -  $(B,  y)]2}  .  (A10) 

For  convenience,  we  introduce  the  notation 

D(A,  6)  i  e{|i|j(z  +  A,  y  ♦  6)  -  iKz,  y )  | 2 }  ,  (All) 

D^(A,  6)  =  e{(x(z  +  A,  y  +  6)  -  x(z»  y))2}  »  (A12) 

D^(  A,  6)  =  e{[$(z  +  A,  y  +  6)  -  4>  ( z ,  y )  ] 2  }  ,  C  A 1 3 ) 

B^(A,  6)  =  E { [ x ( z  +  A,  y  +  6)x(z,  y)  “  e{x(z  +  A,  y  +  6)} 

•  E{x(z,  y)> 1 }  , 

(A  14) 

D  (A,  6)  =  e{(x(z  +  A,  y  +  6)  -  x(z,  y))l<J>(z  +  A,  y  +  6) 

X<t> 

-  $  ( z ,  y)  ] } . 

(A15) 

With  this  notation  and  some  algebraic  manipulation,  (AlO)  becomes 
E{[g(0,  <d;  y)  -  E(g(B,  <d;  y ) ) ) 2 }  =  -D^((D,  0)  -  D^(id,  0)  +  4 

(A16) 


where 


( A 1  7  ) 


l 

> 


f 


°x2  =  e{(x(6,  y)  -  x 1 2 } - 


It  has  been  shown  [6,  12 1  that  from  conservation  of  energy  consid- 

2 

erations,  a  -  _x •  We  also  observe  that 


D( A,  6)  =  D  ,  (A,  6)  +  D  (A,  6). 
9  X 


(A18) 


With  these  substitutions,  we  have 


E{lg(e,  oo;  y)  -  E(g(e,  w;  y))]2}  =  ~D(w,  0)  -  4x*  (A19) 


Finally,  from  (A7 ,  A8 ,  and  A19), 


E{exp[vJ/(0  +  w,  y)  +  i(i*(B,  y)]}  =  exp 


(A20) 


Hence 


oo 

f\  1 

E{  s(x,  y)}  =  |  A  |  J  j  exp  [-  —  D(u>,  0)]  exp[-jku>  •  (x  -  y)] 

-OO 

00 

•  J  J  a(B  +  w)  a(B)  exp[j0(8  +  u)  -  j9(B)]  d2B  d2u>. 

—  OO 

(A21) 

To  compute  the  second  order  statistics,  we  must  develop  an 
expression  for  e|s(xj,  y^)  s(x2>  y j ) } -  From  (9), 
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E{s(v  y2)l  =  I  a  l 4  //  //  //  //  p(b1,  tJ1,  e2,  w2) 

—00  —05  — CO  —00 

2  2 

•  m(bi,  Wj,  yt>  B2,  t»2,  y2)  d  Pj  d  B2  exp[-jka)1  *  ^  -  yj 

2  2 

-  jkw2  •  (x2  -  y2)]  d  a>1  d  0)2  (9) 

where 

PC&1 »  V82»  “  a^Bl  +  “lJ  a^6l^  a^62^  a^62  "  w2^ 

•  exp[je(Bi  +  c^)  -  j e ( 6^ )  +  je(e2)  -  je(e2  -  <*>2)] 

do) 

and 

M(81>  ^ «  S2,  u2,  y2)  =  E{exp[^(61  +  (^ ,  yt )  +  yj 

+  <Kb2,  y2)  +  **(b2  -  u>2,  y2)]}.  (11) 

The  expectation  in  (11)  can  be  computed  by  methods  analogous  to 
those  used  to  compute  the  expectation  in  (A7 ) .  The  calculations 
are  very  tedious,  however,  and  will  not  be  included  here.  Rather, 
we  refer  the  reader  to  Fante  [3]  who  argued  that 
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E|exp[t(Bi  +  U11 ,  ,  y{ )  +  y2 )  +  ^*02  -  u>2>  y2)]} 

=  exp[-  -i-  D^,  O)  -  j  d(w2,  o) 

-  i  D(61  -  e2  +  0.1  +  «2.  yL  -  y2)  -  i  -  &2,  yL  -  y2) 

+  J  d(p1  -  32  +  i»v  yx  -  y2)  +  d(B1  -  B2  +  u2,  yL  -  y2) 

+  2Vei  ~  B2  +  wl»  yl  “  y2^  +  2Bx^l  "  62  +  U2*  yl  "  y2^ 

+  jDx^Bl  '  B2  +  V  yl  "  y2^  '  jDX^Bl  ‘  B2  +  *V  yl  "  y2^- 

(A22) 


Fante  [3]  has  shown  that  the  B  and 

X 

conclude  that 


terms  may  be  neglected. 


We 


m(bj,  ‘»1,  yx,  B2»(02’  y2^  =  expi~  °)  +  d(U2>  °) 

+  d(b1  -  &2  +  <*l  +  u2»  yx  -  y2) 

+  d(Pl  -  P2*  yl  ~  y2 ^ 

-  d(B1  -  P2  +  w1,  yx  -  y2) 

-  d(B1  -  32  +  «2,  yL  -  y2)]}  .  (A23) 
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APPENDIX  B 


In  this  appendix  we  compare  Kummer’s  function  with  the  qua¬ 
dratic  and  5/3-power  approximations.  We  recall  that  Rummer' s  func¬ 
tion  is  the  confluent  hypergeometric  series 


jF1 (a ,  b,  z)  = 


(a)  (a 

(b) (b 


1)  2 
1)  2! 


We  are  particularly  interested  in  the  Rummer's  function  that  ap¬ 
pears  in  (16) 


where  *  5.91/1q  and  1Q  is  the  inner  scale  size  of  the  turbulence 
(typically  1  mm). 

The  quadratic  approximation  of  Rummer's  function  minus  one  is 


!> 


-kV\ 

-H-1 


i- 1) 


v2  2' 
-R  P 
m 


The  term  "quadratic"  approximation  is  derived  from  the  p  term. 
Similarly,  the  5/3-power  approximation  is  given  by 


1  ‘  riii'A.)' 


/kVV6 

I  m  l 

\  4  / 


1 

=  TTTi  i/6) 


IpI  »/ 

m 

(B2) 
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A  comparison  of  Kummer's  function  with  the  quadratic  and  5/3 
power  approximations.  Curve  1  is  Kummer's  function,  curve  2 
is  the  quadratic  approximation,  and  curve  3  is  the  5/3-power 
approximation. 


Rummer's  function  was  computed  numerically  by  summing  the 
first  one  hundred  terms  of  the  series.  This  partial  sum  was  ob¬ 
served  to  converge  over  the  range  of  p  considered.  The  value  of 
Rummer's  function  is  plotted  in  Fig.  10  along  with  the  quadratic 
and  5/3-power  approximations.  The  quadratic  approximation  is  ob¬ 
served  to  diverge  rapidly  from  Rummer's  function  for  p  >  8/Kj^. 
However,  the  5/3-power  approximation  appears  to  be  a  reasonable 
approximation  to  Rummer's  function  for  all  p.  Hence,  we  have  mod- 
eled  jFj(-5/6,  1,  -(r  p  )/A)  by  (23)  throughout  this  paper. 
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APPENDIX  C 


In  this  appendix  we  consider  the  problem  of  select 
coefficients  Cj  and  cQ  in  (51)  so  that 

d(p)  «  c.p  •  p  +  c  |p|  *  w  . 

1  o  c 

For  convenience  we  define 


g(p)  =  c  p  •  p  +  c 
I  o 

We  would  like  the  approximating  function  g(p)  to  be  optimal 
sense.  In  particular,  we  have  considered  four  criteria  for 
ing  the  coefficients  Cj  and  cQ: 

1.  Set 

exp  {-d(p)}jp=u)  =  exp  {-g(p)>jp=aj 

c  c 

and  cQ  =  0. 

2.  Set 


exp  { -d ( p ) }  =  exp  {— gC p) } 

p— (0  p  — CO 


and 


dp  leXp  {_d(p)}1|p=w  =  df  leXP  {_g(p)}1|p=u)  ' 

c  c 


ing  the 

(51) 


(Cl) 

in  some 
select- 
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3.  Minimize  the  functional 


J 


1 


00  oo 

j  exp  (-ad(p)}  dp  -  J  exp  {-ag(p)}  dp 
0  0 


where  a  is  some  positive  constant  (a  >  0). 


4.  Minimize  the  functional 


03 

J  =  JC  [ d ( p)  -  g(p)]  dp 
0 


The  constants  Cj  and  cQ  can  easily  be  solved  for  in 
criteria  listed  above.  The  results  are  summarized  in  Table  5. 


Table  5.  Equations  for  cQ  and  Cj  resulting  from 
each  of  the  four  criteria. 


each 
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From  the  table  we  observe  the  expressions  for  cQ  and  to  be 
reasonably  consistent  for  the  four  criteria  we  have  considered.  We 
have  selected  the  following  expressions  for  c^  and  cq: 


c 


1 


(74) 


c  «>  0.  (75) 

o 
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APPENDIX  D 


In  this  appendix  we  will  describe  in  more  detail  the  numeri¬ 
cal  integration  routine  used  to  integrate  (34).  The  integrand  of 
(34)  is  M' (w,  Aw,  AB,  Ay).  This  function  can  be  loosely  regarded 
as  approximating  a  Gaussian  surface  with  peak  near  AB  =  -Aw/2.  The 
function  rapidly  approaches  zero  for  AB  such  that 

I  A3  +  Aw/2 1  >  w^ .  Hence,  we  should  concentrate  the  numerical  ef¬ 
fort  in  this  region. 

The  numerical  analysis  was  performed  only  for  one-dimensional 

problems.  It  was  desirable  to  use  an  algorithm  that  would  find  the 

peak  described  above  and  automat ical ly  concentrate  the  number  of 

function  evaluations  in  the  vicinity  of  the  peak.  Since  |Aw]  <  w^ 

and  the  width  of  M’(w,  Aw,  AB,  Ay)  in  the  A B  direction  is  about 

w  we  truncated  the  integration  over  the  AB  line  to  the  inverval 
c 

[-1,  +1].  Since  w^_  ~  0.01  to  0.1,  these  bounds  on  AB  are  very 
loose.  To  guarantee  that  the  narrow  peak  of  M'(w,  Aw,  AB,  Ay)  will 
not  be  missed  by  the  routine,  the  interval  [-1,  +1]  was  divided 
into  twenty  subintervals.  An  adaptive  Romberg  [13]  integration 
algorithm  was  applied  to  each  subinterval  and  the  results  were 
summed  to  yield  the  solution.  In  this  manner  the  algorithm  wastes 
little  time  calculating  function  values  of  M' (w,  Aw,  AB,  Ay)  at 
values  of  AB  for  which  the  function  is  small  and  not  rapidly  chang¬ 
ing.  The  computational  effort  is  concentrated  near  the  peak  of 
M'(w,  Aw,  AB,  Ay). 

The  routine  was  found  to  be  robust  but  rather  inefficient. 
The  emphasis  here  was  placed  on  reliable  results  and  not  computa¬ 
tional  efficiency. 
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APPENDIX  E 


In  the  one-dimensional  or  slit  aperture  case,  all  vectors  are 
2  ~ 

colinear  so  cos  8=1,  and  the  model  for  M  (u>,  Aw,  Ay)  can  be 

o 

simplified  to: 


For  the  double  impulse  model,  we  can  extend  the  one-dimen¬ 
sional  model  (El)  into  region  (c)  (section  5)  by  modeling  the  Ay 
dependence  by  the  sum  of  a  constant  and  a  Gaussian.  In  this  case, 
it  is  readily  shown  that 

lira  M  (to,  0,  Ay)  =  --- — —  exp  {-2.92  k^ct  | to |  ^ } 

|h1Ay|?«  {  k  C1ao 

If  we  assume  that  p  =  a  ,  we  obtain  the  model 
o  o 
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1.  Introduction 


The  efforts  reported  here  are  a  continuation  of  the  Speckle  Imaging 
efforts  previously  reported  [1].  Our  goal  was  to  develop  a  Speckle 
Imaging  algorithm  for  non-i soplanatic  conditions.  In  particular,  we 
wanted  a  general  framework  based  on  linear  algebra  in  which  to  examine 
the  many  various  alternatives. 

The  remainder  of  this  report  is  divided  into  nine  sections.  Section 
Two  briefly  reviews  the  non-i soplanatic  speckle  imaging  problem  under 
simplified  conditions.  The  Third  Section  develops  both  a  linear 
algebra  representation  of  the  discretized  speckle  imaging  equations 
and  a  corresponding  graphic  representation.  The  decomposition  of 
rectangular  matrices  is  reviewed  briefly  in  Section  Four.  Section 
Five  develops  the  relationship  that  allows  two-dimensional  images  to 
be  treated  by  the  linear  algebraic  methods  being  developed.  Section 
Six  develops  and  examines  algorithms  for  non-isoplanatic  imaging  based 
on  the  factorization  of  the  second  order  statistics.  The  joint 
estimation  of  the  object  and  the  non-i soplanation  characteri sties  are 
discussed  in  Section  Seven.  Sections  Eight  and  Nine  discuss  the 
effects  of  the  simplifying  assumptions  which  made  the  linear  algebraic 
development  tractable.  Finally,  a  class  of  restoration  algorithms  for 
non-isoplanatic  Speckle  Imaging  is  given  in  Section  Ten. 
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2.  Statement  of  the  Problem 


Our  aim  In  this  work  is  to  restore  telescopic  images,  namely  to  find 
the  original  image,  the  object,  from  images  corrupted  by  the  turbulent 
atmosphere.  The  work  is  based  on  an  analysis  by  Sherman  [1].  In 
particular  we  wish  to  consider  the  removal  of  what  is  called  the 
non-isoplanatic  effect  through  use  of  second  order  statistics  namely 
the  average  of  many  identical  images  (save  atmospheric  effects)  of 
some  celestial  object.  The  above  paper  provides  the  analyses  and  we 
take  as  a  starting  point  its  equation  38  which  we  want  to  solve  for 
the  original  object.  It  behooves  us  here  to  try  to  give  some 
intuitive  interpretation  to  this  equation  or  rather  to  the 


corresponding  relationships  of  the  Fourier  transformed  variables. 
Thus,  if  we  denote  by  a  tilda  ('*T  a  Fourier  transformed  variable,  the 
corresponding  equation  which  we  call  (38)  is  obtained  from  (38)  by 

/v/ 

eliminating  the  integrals  with  the  inverse  Fourier  kernel.  Thus  (38) 
is 


o-' 


—  wit. 


In  (38),  0  is  the  transform  of  the  object,  while  I  is  the  transform  of 
the  image.  To  understand  the  mechanism  let  us  substitute  for  the 
expectation  of  (38) 


namely  the  left  hand  side  of  (38)  is  approximated  by  the  sample, 
second  order  moment  of  Q  images.  Because  the  effect  of  the  turbulence 
could  be  considered  as  adding  a  random  phase  component  in  the  Fourier 
plane,  (1)  would  generally  represent  summation  of  complex  numbers  with 
an  added  random  phase  which  would  cause  partial  to  full  cancellation 
except  in  the  case  where 
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Then  the  product  is  always  positive  with  the  sum  adding  up  because 

ft»i) 

ho l Os  for  an  original  real  (physical)  image.  We  thus  could  say  (proof 
in  [1]  )  that  the  atmospheric  effect  on  the  expection  of  (1)  for 
)m,l  and /^i/  >  u>c  is  to  multiply  what  would  have  been  otherwise 

obtained  by  M  (u.>,4&>),  where  Lo  and  OL)  are  given  (as  defined  in  [1] 

)  by: 


V 


Co, 


(3) 


(V,  * 

x. 


•  LO 
J 


(4) 


The  U),Ab)  plane  (taking  for  illustration  pusposes  only  one  dimension 
of  each  of  it’andU.^)  is  rotated  with  respect  to  the  U),“> z  plane  by 
45°  (and  actually  also  reduced  by  a  factor  of  Y1T ).  We  show  in 
Figure  1.  is,  except  near  the  origin,  a  function  mainly  of  and 
only  weakly  of  L) .  Thus,  on  the  line  is  maximum  and  it 

tapers  off  as  we  move  away  from  this  line.  It  can  be  approximated  by 
a  Gaussian  curve.  We  also  have  a  factor,  F^  ,  which  measures  overlap 
of  the  aperture  with  four  shifts.  What  it  signifies  is  that  the 
Fourier  plane  is  truncated  by  the  aperture  of  the  telescope  and  that 
the  overlap  of  it  appears  in  the  correlation.  F^  is  given  by 


where  a(*)  is  one  inside  the  aperture  and  zero  outside.  We  may  note 
that,  because  we  are  interested  only  in  the  region  around  6/=-^ where 
HjfO  and  that  a(')=a*t'),  (5)  reduces  to  the  two  shift  overlap 
function  for  the  region  of  interest. 


represents  in  (38)  the  non-isoplanatic  effect.  This  is  the  effect 
we  want  to  remove,  or  more  precisely,  the  effect  despite  which  we  want 
to  achieve  reconstruction  of  the  object  from  Fourier  plane  data.  The 
impluse  response  of  the  imaging  system  and  atmosphere  combined  is,  in 
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► 

I 


system  theory  terminology,  shift  variant  under  non-i soplanatic 
conditions.  The  interpretation  of  (38)  means  that  the  effect  of  Mjis 
that  the  value  of  the  expected  cross  correlation  at  a  point  (  ^,^2) 
is  given  apart  from  the  effect  by  smoothing  along  a  segment  of 
constant  with  as  a  weighting  function.  Figure  1  illustrates 
this  smoothing  graphically.  The  weighting  function,  Mj,  is 
approximately  a  Gaussian  function  with  a  width  of  k.^. 
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3.  Algebraic  Formulation  of  Discretized  Equations 


We  rewrite  (3S)  here  as 


(6) 


where  the  (~!  notation  has  been  dropped  and  &  is  taken  only  in  the 
range  of  support  of  and  F^.  Let  us  assume  for  awhile  that  LOl  and 
LOx  are  one  dimensional  (rather  than  two).  Let  us  assume  also  that  we 
digitize  (6)  on  a  uniform  grid.  Then  we  have  a  discrete  version  of 
(6), 


(mumx) 


y  0(nnl  -J)  Of  d)  AJ(J ) 


(7) 


0(m)  is  a  (column)  vector  and  thus  we  can  consider  O(m-k)  to  be 
another  vector  generated  from  the  first  by  a  shift  of  elements.  Thus 
all  possible  values  of  O(m-k)  for  L<m<H,  -R<k<R  could  be  represented 
by  a  Toeplitz  matrix  whose  elements  are  O(m-k). 


7 


o  = 


m 


Of  L-tf\) 


0(l~ X) 


(8) 


1  O(nj-J) 


O(H-n) 


where  the  left  subscript  T  means  that  the  Toeplitz  matrix  r0  is 
obtained  from  the  vector  0  by  appropriately  shifting  the  elements. 

Because  Mj  is  similar  to  the  Gaussian  function,  we  select  R  according 
to  the  precision  required.  Similarly  we  can  define  yyO  as  the 
generation  of  a  Hunkel  matrix  whose  elements  are  0(m+k)  from  the 
vector  0.  Note  that  while  the  elements  of  the  Toeplitz  matrix  arc1  the 
same  as  diagonals,  they  are  the  same  on  anti  diagonals  in  Hunkel 


-6- 


matrices.  Thus  (7)  can  be  written  as 


-  TonHo  (9) 

L_  J  '  " 

where  superscript  T  indicates  the  transpose  and  M  is  a  diagonal  square 
matrix  whose  elements  are 

•=  Mi 

Figure  2  graphically  represents  (9).  Any  element  ofo^,  located  say 
at  m,n  is  obtained  by  a  scalar  or  inner  product  of  row  m  from  T0  times 
column  n  of  ^0  (or  row  n  of  ^0),  weighted  by  the  diagonal  elements  of 
M. 


We  should  also  note  that  a  Hunkel  matrix  could  be  obtained  from  a 
Toeplitz  matrix  and  visa  versa  by  multiplication  by  an  '‘antidiagonaT’ 
unit  matrix 


so  that  (9)  would  become 

=  jONJjO  U2) 

We  note  also  that  in  (9)  or  (12),  M  could  be  absorbed  into  r0  and 
being  diagonal  M  multiplies  each  column  of  r0  by  a  constant  mjj .  It 
of  course  could  be  also  absorbed  into  the  right  as  M=MT,  or  we  could 
split  M,  say  symmetrically,  into  M=7rM]tM.  Also  M  commutes  with  J  to 
produce  for  example, 

z  7Offi  J~ llMjO  (13) 
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4.  Decomposition  of  Rectangular  Matrices 


We  wish  first  to  review  and  present  in  a  form  to  be  used  later  some 
results  from  linear  algebra.  Consider  any  pxq  rectangular  matrix 
A=[a(i,j)]  of  rank  r.  This  implies  that  we  can  find  r  independent 
columns  of  this  matrix  and  express  all  other  columns  as  linear 
combinations  of  these  columns.  Thus  we  can  write 


A  -  /)*  Al  . 


where  A^is  the  pxr  matrix  generated  by  the  above  columns  and  AL  is  a 
qxr  matrix  whose  columns  are  the  coefficient  of  the  linear  combination 
that  generate  the  columns  of  A.  We  could  alternatively  start  by 
saying  that  any  row  of  A  is  expressed  as  a  linear  combination  of  r 
rows.  Thus,  it  is  said  that  A«  and  are  respectively  composed  of 

^  T 

vectors  which  span  and  are  a  basis  of  the  columns  space  of  A  and  A  . 
This  decomposition  is  called  "full  ran',  decomposition  or  "full  rank 
factorization".  The  name  implies  the  full  column  rank  of  A^  and  row 
rank  of  A 7".  This  is  shown  graphically  in  Figure  3.  We  should  note 
that  the  decomposition  is  not  unique,  because  all  possible  A^  span  the 
same  column  space  of  A  and  all  At  span  the  row  same  space  of  A.  One 
A^  is  related  to  another  by  a  non-singular  transformation,  and  we 
actually  can  also  write  a  more  general  decomposition 

A-zeeJ,  ns) 


where  C  is  a  non-singular  rxr  matrix  (Figure  3c).  We  may  note  that  we 
could  use  a  decomposition  with  a  number  of  columns  in  B*  and  rows  in 
B^that  is  larger  than  the  rank.  This  is  more  general,  but  we  will 
lose  some  properties. 


While  it  may  seem  trivial,  it  is  worthwhile  to  examine  Figure  3  in 
detail.  Considering  Figure  3b,  we  see  that  each  element  of d  ,  say 
is  obtained  by  a  scalar  or  inner  product  of  the  corresponding  row 

Q 


W  1  11  11  "I— 


r 

from  and  of  the  corresponding  column  of  kL.  Looking  at  Figure  3C, 
we  see  that  we  can  absorb  C  into  Bp.  This  generates  a  new  in  which 
each  column  is  a  linear  combination  of  those  of  the  old  B^.  One  could 
of  course  absorb  it  in  and  then  the  operation  would  be  on  the  rows. 
C  could  also  be  split  into  a  product  C;  and  then  C,  would  operate  on 
the  columns  of  B^  and  on  rows  of  B^  .  In  the  particularly  useful 
case  when  C  is  diagonal,  what  happens  is  that  each  column  of  A^is 
multiplied  by  the  corresponding  diagonal  element  which  is  a  constant. 

We  should  also  note  that  the  matrix  A  can  be  really  only  a  submatrix 
of  some  larger  matrix  that  has  been  partitioned.  If  we  consider 
Figure  4  or  the  corresponding  equations 

We  note  that  we  cannot  take  all  i  or  j  but  only  some  of  them.  Each 
element  in  sJ  is  obtained  from  a  row  B^  and  a  column  of  B^.  Thus  an 
independent  subset  of  equations  can  be  obtained  if  we  select  a  set  of 
rows  and  columns  from  B^  ,bJ"  as  shown  in  Figure  4.  In  fact  the 
selected  rows  or  columns  need  not  be  contiguous.  One  should  be 
careful  to  note  that  this  selection  may  change  the  nature  of  the 
equations  represented  by  the  matrix  product.  For  example,  B^,  and 
B^  may  contain  the  same  variable  which  would  make  the  equations 
quadratic.  Proper  selection  may  change  this  nature.  This  possibility 
of  selection  is  important  in  particular  cases  both  because  it  allows 
reduction  in  the  number  of  equations  (at  the  cost  of  redundancy  or 
noise  reduction)  and  change  in  nature  of  relationship. 

A  particular  form  of  full  rank  decomposition  that  has  become  very 
popular  in  recent  literature  is  the  singular  value  decomposition 
(SVD).  The  SVD  is  a  full  rank  decomposition  with  two  additional 
requirements.  First,  the  diagonal  matrix  C  in  (15)  of  positive 
elements  (usally  called  2_)  is  arranged  by  magnitude  with  the  largest 
one  first.  Second,  B  and  B  (usually  called  U  and  V  respectively) 

A 
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are  composed  of  orthogonal  vectors  of  unit  length.  The  SVD  is  usually 
written  with  the  notation 

fy  -  U  £  V  or  /j  --  2  Li -7,1:  (17) 

c  -  / 

if  A  is  square  and  Hermitian  then  U  is  the  conjugate  of  V. 

In  some  presentations  in  the  literature  U  and  V  are  taken  to  be  square 
rather  than  of  a  number  of  columns  equal  to  r,  the  rank  of  A.  This  is 
done  to  be  able  to  make  U,V  conform  to  the  definitions  of  orthogonal 
matrices.  Because  the  u£ ,  for  i  such  that  OJ  =0  do  not  affect  A, 
the  two  definitions  coincide.  Because  the  SVD  is  a  full  rank 
decomposition  (factorization) ,  we  may  conclude  that  one  can  represent 
that  matrix  A  by  either  the  SVD  representation  or  the  full  rank 
decomposition  of  (15)  and  that  the  transformation  from  U  to  B^  (or  V 
to  BL)  is  a  nonsingular  transformation  which  is  recoverable  if  we  know 
an  r  row  section  of  U  and  B^ . 
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5.U  Representation  of  Two-Dimensional  Data 

We  previously  developed  the  analysis  on  the  simplified  assumption  that 
our  object  and  image  were  one  dimensional  and  this  led  us  to  a  matrix 
representation  for  the  basic  non-i sopl anati c  relation  for  second  order 
statistics.  Because  our  objects  are  two  dimensional,  their  Fourier 
transforms  are  also  two  dimensional.  We  can  write  the  sampled  version 
as  a  function  of  two  discrete  variables,  n,m  or 

0  ^  0(r*),n]  (i9) 

and  two  dimensional  version  of  (7),  the  basic  equation,  is 


and  our  functions  now  have  four  variables  rather  than  two. 


Because,  we  are  used  to  representations  of  functions  of  two  variables 
and  they  are  easier  to  represent  in  the  plane  (paper)  graphically,  we 
would  like  to  map  the  four  dimensional  Euclidian  space  into  a  two 
dimensional  one.  This  can  be  done  in  several  ways.  We  start  with  the 
way  we  look  at  the  original  problem.  Figures  5a  and  5b.  We  have  a  2-d 
plane,  and  each  point  in  it  is  the  origin  of  a  perpendicular  2-d 
plane,  but  the  axes  can  be  selected  in  several  ways.  We  deal  not  with 
continuous  variables  but  with  a  finite  set  of  discrete  variables. 

These  can  be  mapped  from  a  plane  to  a  line  by  taking  a  column  at  a 
time  and  cascading  them.  In  order  to  be  consistent  with  our  previous 
representations,  we  will  use  the  mapping  ofc^  illustrated  in  Figure 
6.  The  index  space  will  map  into  a  column  by  taking  a  concatenation 
of  vectors  running  over  n  for  fixed  m.  This  mapping  is  the  common 
raster  scan  and  keeps  (20)  in  its  familiar  form, 

J  --tOMhOT  <«, 
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where  now  the  matrices^,  _Q,  .0  become  block  matrices.  Thus,  _0 

in  T 

could  now  be  considered  as  block-block  Toeplitz.  This  means  that  in 

the  m,k  block  location  of  matrix  T0  we  have  a  Toeplitz  block  with 

indices  n,l  and  the  block  themselves  have  a  Toeplitz  structure. 

Similarly,  ^0  is  block-block  Hunkel  and  has  a  corresponding  block 

sturcture.  M  is  now  again  a  diagonal  matrix,  but  the  values  of  the 

elements  are  Mf-vV+l2  )  and  do  not  have  the  bell  shape  any  more. 

Figure  7  shows  the  mapping  for  theT0  matrix  for  an  arbitrary  vicinity 

of  k,l;  1  ranging  from  -1  to  +  1  (assuming  R=l).  This  makes  the 

width  of  the  T0  matrix  9=Niwhere  N=(2R+1).  If  the  m  ,n  array  has 

dimensions*  uxv,  then  T0  is  a  block  matrix  of  uxN  blocks  and  each 

'  % 
block  is  vxN.  The  total  matrix  dimension  is  uvxN  . 


*Please  note  the 
and  in  a  space. 


different  meaning  of  the  word  dimension  in  a  matrix 
The  appropriate  meaning  will  be  clear  from  context. 
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Figure  7.  Mappings  of  two-dimensional  indices  to 


6.  Algorithms  for  Factorization 


We  will  introduce  here  an  algorithm  for  finding  0  given^and  Mj.  We 
show  in  Section  7  that  if  Mj  is  not  known  it  can  be  found  as  part  of 
the  procedure.  For  this  purpose  let  us  consider  first  the  Newton 
Raphson  procedure  for  finding  a  square  root  or  finding  X  given 


XZ  =  A 


Let  be  the  1th  iterate  obtained  using  this  Newton  Raphson  formula 

4,  <23> 


X,.,  =  i  (ti '%) 


We  purposely  separated  (23)  and  (24).  While  (24)  seems  more  compact, 
it  may  be  advantageous  numerically  to  use  (23). 

These  formulas  are  known  to  converge  very  strongly-  namely 
quadratically.  A  simple  substitution  shows  that 


—  'ty-f/  X  J-( 

x  v  x  l 


X  "  — -  ^  l  x  / 

So  if  we  have  a  fair  approximation,  say  .25X,  then  after  the  next 

iteration  we  have  X,  -  3%X. 

Jn 


Analogously  we  start  with 


J  ---pMHoT. 


We  take  the  left  inverse  of  T0  and  obtain 


oLj -  (o' o)'dJ --  Nh o':W° 
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-19- 


where  again  J  is  the  "anti diagonal  unit"  matrix. 


We  now  assume  that  ^  is  the  1th  approximation  to  r0  and  define 

U*l  =  (AlDi]70!^  (281 

We  also  define  T~J 0  as  the  (n+p-1)  vector  obtained  from  the  nxp  matrix 
T0  by  averaging  all  the  elements  on  the  diagonals.  jO  may  not  be 
exactly  Toeplitz,  or  as  is  the  case  here  the  columns  may  be  multiplied 
by  the  factor  .  We  define  h-i  0  similarly.  Now  we  let 


vhi  -zMu 


A,„ 


where  the  sum  is  only  over  the  terms  in H  corresponding  to  the 
elements  averaged.  This  is  the  complete  range  of  k  for  diagonals 
which  are  not  foreshortened. 


Now  let 


a  A,  -  f  Tl^"~  °A 


A  ATl 


tOj„  A  -[  Vtni 


the  same  comments  about  numerical  advantage  apply  to  (31)  and  (32)  as 
those  of  (23)  or  (24).  Straight  application  of  this  algorithm  assumes 
Jj  (or  rather  the  selected  rows  and  columns  of gj  )  is  square. 


A  variant  of  this  procedure  which  was  used  in  [13  is  to  normalize 
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(multiply  by  constant)  in  such  a  way  that  the  (quadratic)  norm 
implied  from  (26)  will  hold.  This  converges  for  the  rank  one  case  in 
one  iteration.  Consider  a  matrix  C  of  rank  one 

C  ~  X  X  T  (33) 

We  let 

x'~  Cl  =  xXxt ■ =  X/x\)  (34) 

and 

\  -  X  1/Co  (35) 


where  the  scaling  constant  c0  is  defined  by 

r  4  C  l 

and  the  matrix  dot  product  is  defined  by 

=  A‘j  4 


Then  c0  can  be  reduced  to 


Cc 


Xe 


7 


(36) 


(37) 


(38) 


Substituting  into  (35)  shows  that  (34)  and  (35)  converge  in  one 
iteration,  when  C  is  rank  one, 


T  /  /7V  f  1 — 

X(xk)(xx«)  - 


X 


(39) 


If  there  Is  more  than  one  non-zero  eigenvalue,  then  convergence  will 
not  happen  in  one  step.  The  error  after  the  1th  iteration  decreases 
as  the  ratio  of  the  second  largest  eigenvalue  to  the  largest 
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eigenvalue  raised  to  the  1th  power.  This  method  is  an  adaptation  of 
the  power  method  for  computing  eigenvectors. 
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1 


7.  Joint  Estimation  of 


We  assumed  that  Mj  is  given  to  a  first  approximation  by 

M3( •»')--  «-  ifc)  (40' 

Actually  it  turns  out  that  both  ^'c  and  a.  are  weak  functions  of  cu  ,  as 
shown  in  Figure  8.  As  seen  from  the  figure,  we  can  divide  the 
behavior  of  Mj into  three  regions.  In  region  I,  low  frequencies,  it 
is  essentially  an  impulse.  Region  II  is  a  transition  region  in  which 
it  widens  to  a  guassian  funcion  and  becomes  fairly  independent  of  in 
region  III.  Because  the  low  frequency  behavior  of  0  is  known,  our 
main  interest  is  in  region  III  and  we  have  modelled  the  behavior  of 
Mj  in  another  portion  of  this  effort. 


If  we  do  not  know  Mj,  we  can  find  it  as  part  of  our  computation  to 
find  0.  As  we  have  seen,  we  obtain  a  Toeplitz  matrix  and  each  column 
of  it  is  multipled  by  a  sampled  value  of  Mj  which  we  called  M \j^.  Thus 
at  one  of  the  computational  stages  we  may  obtain  (noting  _^=M ^ ) 

the  matrix 


(41) 
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Because,  apart  from  ,  elements  on  each  diagonal  are  the  same,  it 
seems  that  we  can  take  any  diagonal  and  find  M^up  to  a  constant. 

(Mm  is  a  positive  number).  However,  the  entries  may  be  corrupted  by 
noise  or  numerical  inaccuracy. 

If  we  look  at  (41),  we  see  that  not  only  elements  on  the  diagonals 
differ  only  by  their  multiplier  but  also  the  encircled  elements  as 
a  group.  At  first  glance  it  would  seem  advantagous  to  sum  all  the 
elements  in  encircled  groups  and  normalize  the  results  so  that  =1 . 
The  Oj  are  however  complex  numbers  and  their  sum  may  cancel.  It  would 
behoove  us  thus  to  take  the  sum  of  the  magnitudes  in  our  example  or  to 
evaluate 


J, 


The  number  of  elements  in  the  summation  (a+b+1)  is  chosen  according  to 
considerations  numerical  of  accuracy,  amount  of  computation,  and 
signal-to-noise  ratio. 
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8.  Effects  of  Truncating  of 

We  have  stated  that  we  take  R  such  that  Because  Mj  only  goes 

to  zero  assymptotical ly ,  though  very  strongly,  we  would  like  to  find 
the  effect  of  truncating  the  "tails'1  of  Mj.  A  smaller  R  is  desirable, 
because  it  reduces  the  number  of  equations  (proportional  to  R  )  and 
improves  the  condition  nutroer. 


Figure  9  illustrates  the  product 

C-  A  h  E‘ 


(44) 


where  A  is  mxq,  diagonal  M  is  qxq,  and  B  is  nxq.  We  want  to  find  the 
change  in  the  elements  of  C  if  q  is  taken  to  be  one  less.  We  denote 


by 


r- 


Ai  c. 


(45) 


the  change  in  the  elements  of  C  generated  by  truncating  the  qth 
element  of  Hj.  The  term  in  brackets  is  a  rank  one  matrix  or  dyadic. 
Because  A,B  are  respectively  Toeplitz  and  Hunkel  matrices,  all  their 
columns  have  about  the  same  norm.  Thus,  the  individuals  terms  are 
proportonal  to  and  the  effect  of  truncation  of  the  sum  at  a 
certain  q  is  proportional  to  the  "area"  under  the  tails  of  M  which  is 
approximately  Gaussian  and  tapers  off  rapidly. 
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9.  Effects  of  the  M,  Schur  Product 


As  indicated  before  MjF^  multiply  each  term  of  *=/.  This  is  an 
operation  known  algebraically  as  "Schur  Product"  and  which  does  not 
work  well  together  with  the  ordinary  matrix  multiplication.  One  could 
divided  by  MjFv,  but  this  accentuates  the  noise  and  is  not 
desirable.  The  basic  non-isoplanatic  equation  (38)  can  be  written  as 


^  v  c  Fi  c  (-fi  ^  )■  la?) 


One  way  to  overcome  this  difficulty  is  to  consider  only  one  row  at  a 
time  of  MjFy  .  Then  we  could  take  that  one  row  of  MjF^  and  write  it  as 
a  diagonal  matrix  for  one,  particular^,.  This  is  depicted  in  Figure 
10  which  is  similar  to  our  previous  representations,  with  added  and 
denoting  the  multiplication  of  each  column  of^O  by  the  corresponding 
(diagonal)  element  of  M ^  , 


OMhOt^ 


1 


(48) 


We  have  absorbed  Fy  into  M,  and  we  can  solve  for  one  row  of^O, 

ra=j<o(vHOTfytf* 

where  -R  indicates  a  right  inverse. 


Thus  we  can  find  one  row  ofr0  at  a  time.  Because  r0  is  a  Toeplitz 
matrix,  we  do  not  need  all  rows  to  find  the  vector  0  but  in  order  to 
get  more  redundancy  we  would  compute  one  row  out  of  about  every  ui'd 2 
rows.  The  amount  of  work  is  large,  because  we  compute  every  right 
inverse  anew.  Although  this  difficulty  could  be  resolved  if  M*  were 
separable,  is  not  separable. 
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10.  Restoration  Algorithm 


We  are  now  ready  to  write  a  restoration  algorithm  for  the 
non-isoplanatic  case  in  simple  algebraic  notation.  The  Speckle 
Imaging  equations  for  non-isoplanatic  conditions  can  be  written  as 
functions  of  two  dimensional  variables  ,  ay  and  u/ , 

Alfa.n'i) 

'  lo/'K< 

Again,  and  have  been  combined  into  M.^  We  will  assume  that  the 
sample  second-order  statistics  denoted  by  c=P  (k’(  ,&' )  have  been  corected 
for  the  bias  terms  introduced  by  measurement  noise. 

A 

Fortunately , (w, ,Aa)  does  not  have  to  be  computed  for  all  a’;  and  ay 
The  second-order  statistics  have  non-negligible  values  and  hence 
information  about  the  object  only  for  Thus,  we  only 

crosscorrelate  those  portions  of  the  image  transforms  which  are 
mutually  coherent.  The  size  of  the  mutually  coherent  area  is 
described  by  the  atmospheric  cutoff  frequency,^  .  The 
non-isoplanatic  effect  introduces  a  smoothing  of  the  second-order 
statistics  as  indicated  by  the  summation  over^/in  (50).  This  allows 
a  further  reduction  in  the  number  of  ^ltlVx  points  at  which  the 
second-order  statistics  must  be  computed  because  the  redundant  values 
do  not  contribute  substantially  to  the  solution. 


As  indicated  previously,  we  can  solve  for  the  object  by  solving 
subsets  of  equations  from  the  set  of  equations  given  by  (50).  For 
example,  select  all  equations  involving  a  fixed  UJt  and  write  them  as  a 
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Because  the  isoplanatic  patch  is  greater  than  the  spread  response, 
less  than&£.  This  results  in  (51)  having  more  equations  than 
unknowns  because  there  are  more  values  of  Iuj^I^u^  than  This 

formulation  is  illustrated /in  Figure  10  as  the  selection  of  a 
particular  row  of  T0  and  ck.  The  remaining  question  of  the  linear 
independence  of  the  equations  was  previously  addressed  in  [1]. 

Because  the  subset  of  equations  is  based  on  a  fix  _  ,  they  are 
independent.  The  equations  can  be  written  in  vector  form  as 


Lo> 


I(v>.)  -  & )  0(^<) 


(54) 


The  solution  can  be  written  as 


However,  the  inverse  would  not  be  computed  because  it  depends  on 
6tJ(  and  the  current  estimate  and  can  be  used  only  once.  Directly 
solving  (54)  using  a  linear  equation  solution  method  would  be  more 
computationally  efficient.  This  yields  an  estimate  of  the  object  over 
a  region  of  radius  ooc  about  W  . 

To  solve  for  the  object  over  the  entire  Fourier  transform  domain,  we 
choose  a  sequence  of  at  which  we  solve  (54).  The  sequence  of 
U>,' s  should  start  at  the  origin  and  "spiral"  out  to  the  diffraction 
limit.  The  separations  between  then's  chosen  should  be  less  than 
twice  Dc  to  allow  the  regions  of  the  estimates  to  overlap  and  cover 
the  entire  frequency  domain,  and  should  be  greater  than  D*  to  avoid 
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processing  redundant  information. 

The  steps  of  this  restoration  algorithm  for  non-isoplanatic  Speckle 
Imaging  can  be  summarized  as  follows: 

1.  Initialize  the  solution  near^,=0  using  the  first-order 
statistics  (average  image  transform). 

2.  Compute  the  coefficients  of  the  linearized  equations,  C (&',), 
from  the  current  estimate  and  the  models  for  atmospheric 
disturbance  and  measurement  noise. 

3.  Crosscorrelate  the  second-order  statistics  with  C (u'f)  and 
solve  (54)  for  the  object. 

4.  Select  a  new  uj  and  repeat  starting  at  step  2  until  the 
diffraction  limit  is  reached. 

5.  Reset <V,  to  zero  and  repeat  starting  at  step  2  until  a 
convergence  criterion  is  satisfied  or  for  a  predetermined 
number  of  iterations. 
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